clear; clc; close all;
% Extract the load demand from excel file
Load = xlsread('Book1.xlsx', 1);
% get the number of hours for optimazation operation
nHours = numel(Load);
Time = (1:nHours);
%Enter the Number of Generating Units in the System
num_of_gen = 6;
plant = cell(1, num_of_gen);
%Enter Generator Data from Excel File
%Create Matrix to store data
gen_data = zeros(num_of_gen, 9 );
gen_data = xlsread('Book2.xlsx', 1);
%Pull apart generator properties
%Get Start Up Cost data
StartupCost = (gen_data(:, 9));
%Calculates F(P_i,min)/Operating Cost
OperatingCost = zeros(1, num_of_gen);
for num = 1:1:num_of_gen
%Operating Cost = a.(P_min)^2 + b.(P_min) + c
OperatingCost(1, num) = gen_data(num, 2) * (gen_data(num, 6) ).^2 + (gen_data(num, 3) * gen_data(num, 6)) + gen_data(num, 4);
end
%Get the Minimum up and down times data
MinimumUpTime = (gen_data(:, 7))';
MinimumDownTime = (gen_data(:, 8))';
%==========PIECEWISE LINEARIZATION OF COST FUNCTIONS==========
%Find Segment Width
%Specify the number of Segments
%Increasing the number of segments improves the results
K = 4;
W = zeros(num_of_gen, 1);
for num = 1:1:num_of_gen
%W = (P_max - P_min)/K
W(num, 1) = (gen_data(num, 5) - gen_data(num, 6)) / K;
end
%Get Maximum generation level for each segment of linearized cost curve
MaxGenerationLevel = zeros(1, (num_of_gen*K));
for num = 1:1:num_of_gen
for k = 1:1:K
MaxGenerationLevel(1, ((num - 1)*K)+k) = W(num, 1);
end
end
%Set Minimum generation level to 0
MinGenerationLevel = zeros(1, (num_of_gen*K));
%Get Maximum generation level for the actual units
MinPowerGenLevel = gen_data(:, 5)';
%Find end point for each Power Segment
Power_Segment = zeros(num_of_gen, (K + 1));
for num = 1:1:num_of_gen
for k = 1:1:(K + 1)
if k == 1
Power_Segment(num, k) = gen_data(num, 6);
elseif k == (K + 1)
Power_Segment(num, k) = gen_data(num, 5);
else
Power_Segment(num, k) = gen_data(num, 6) + ((k - 1) * W(num, 1));
end
end
end
%Find Respective Cost for each end point of each Power Segment
Cost = zeros(num_of_gen, (K + 1));
for num = 1:1:num_of_gen
for k = 1:1:(K + 1)
Cost(num, k) = (gen_data(num, 2) * Power_Segment(num, k) * Power_Segment(num, k)) + (gen_data(num, 3) * Power_Segment(num, k)) + gen_data(num, 4);
end
end
%Find Respective Slope for each Power Segment
Slope = zeros(num_of_gen, K);
for num = 1:1:num_of_gen
for k = 1:1:K
Slope(num, k) = (Cost(num, (k+1)) - Cost(num, k)) / W(num, 1);
end
end
%===================END of LINEARIZATION======================
%======FORMULATION FOR USE OF 'intlinprog' solver FOR OPTIMIZATION====
%Build Objective Function Matrix
f = zeros((num_of_gen * K), 1);
%Create sub matrices to easily formulate f
sub_f = cell(num_of_gen, 1);
for num = 1:1:num_of_gen
for k = 1:1:K
sub_f{num, 1}(k, 1) = Slope(num, k);
end
end
%Convert cell to matrix to formulate f
f = cell2mat(sub_f);
%This code considers each segment of the linearized cost curve as a
%separate plant
plant = cell(1, (num_of_gen*K));
for num = 1:1:num_of_gen
for k = 1:1:K
plant{1, ((num - 1)*K)+k} = ['P' '_' num2str(num,'%d') num2str(k,'%d')];
end
end
%The actual number of plant\generators in the system
actual_plant = cell(1, num_of_gen);
for num = 1:1:num_of_gen
actual_plant{1, num} = ['P' '_' num2str(num,'%d')];
end
%Create maxGenConst, minGenConst and minPowConst
nSlots = nHours*num_of_gen;
idxHr2ToEnd=2:nHours;
maxGenConst = repmat(MaxGenerationLevel,nHours,1);
minGenConst = repmat(MinGenerationLevel,nHours,1);
minPowConst = repmat(MinPowerGenLevel,nHours,1);
%Create an optimization problem
powerprob = optimproblem;
%Create Optimization Variables
%Amount of power generated in an hour by a plant
power = optimvar('power',nHours,plant,'LowerBound',0,'UpperBound',maxGenConst);
%Indicator if plant is operating during an hour
isOn = optimvar('isOn',nHours,actual_plant,'Type','integer','LowerBound',0,'UpperBound',1);
%Indicator if plant is starting up during an hour
startup = optimvar('startup',nHours,actual_plant,'Type','integer','LowerBound',0,'UpperBound',1);
%maxGenConst is a 24x(num_of_gen*K) matrix
%isOn is a 24x(num_of_gen) matrix
%Need to do piece-wsie multiplication later on the update LB and UB
%The following converts maxGenConst to a 24x(num_of_gen) matrix
%And then transforms it back to a 24x(num_of_gen*K) optimization matrix
sub_maxGenConst = zeros(24, num_of_gen);
for num = 1:1:num_of_gen
sub_maxGenConst(:, num) = maxGenConst(:, (num * K));
end
new_maxGenConst = sub_maxGenConst.*isOn;
sub_new_minGenConst = zeros(24, (num_of_gen*K));
%Initialise sub_new_minGenConst
sub_new_minGenConst = repmat(minGenConst(:, 1), 1, K);
for num = 2:1:num_of_gen
sub_new_minGenConst = [sub_new_minGenConst repmat(minGenConst(:, num), 1, K)];
end
%Costs
powerCost = sum(power*f,1);
isOnCost = sum(isOn*OperatingCost',1);
startupCost = sum(startup*StartupCost,1);
%Set objective
%Minimize the following:
powerprob.Objective = powerCost + isOnCost + startupCost;
%Set up Constraints
%Load Demand Constraint i.e. the total power generated must meet load demand
%Power generation over all plants meets hourly demand
powerprob.Constraints.isDemandMet = sum(power,2) >= Load - sum((minPowConst.*isOn),2);
%enforce startup=1 when moving from off to on
% no need to enforce startup=0 at other times since minimizing cost forces it
powerprob.Constraints.startupConst = -isOn(idxHr2ToEnd-1,:) + isOn(idxHr2ToEnd,:) - startup(idxHr2ToEnd,:) <= 0;
%Min uptime Constraint
powerprob.Constraints.minUptimeConst = optimconstr(nHours,actual_plant);
for jj = 1:num_of_gen
for kk = 1:nHours % based on possible startups; no penalty at end for running over
if kk > nHours-MinimumUpTime(jj)
sumidx = kk:nHours;
else
sumidx = kk:kk+MinimumUpTime(jj)-1;
end
powerprob.Constraints.minUptimeConst(kk,jj) = ...
startup(kk,jj) - sum(isOn(sumidx,jj)/length(sumidx)) <= 0;
end
end
%Min downtime Constraint
powerprob.Constraints.minDowntimeConst = optimconstr(nHours,actual_plant);
for jj = 1:num_of_gen
for kk = 2:nHours % based on possible startups; no penalty at beginning
if kk <= MinimumDownTime(jj)
sumidx = 1:kk-1;
else
sumidx = kk-MinimumDownTime(jj):kk-1;
end
powerprob.Constraints.minDowntimeConst(kk,jj) = ...
startup(kk,jj) + sum(isOn(sumidx,jj)/length(sumidx)) <= 1;
end
end
%==================END of FORMULATION=====================
% options for the optimization algorithm, here we set the max time it can run for
options = optimoptions('intlinprog','MaxTime',100000);
% call the optimization solver to find the best solution
[sol,TotalCost,exitflag,output] = solve(powerprob,options);
%Calculate the Total Power Generated
Total_Power = zeros(24, num_of_gen);
% %Create a sub_Total Power to contain P_min based on isOn
sub_Total_Power = minPowConst .* sol.isOn
%Populate Total_Power
%Initialise Total_Power with sub_Total_Power
Total_Power = sub_Total_Power;
for h = 1:1:24
for num = 1:1:num_of_gen
for k = 1:1:K
Total_Power(h, num) = Total_Power(h, num) + sol.power(h, (((num - 1) * K) + k));
end
end
end
%Find the Operation cost for each unit for each hour
sub_Hourly_Cost = zeros(24, num_of_gen);
for h = 1:1:24
fo
