Optimization Problems in MATLAB

I was looking at solving the Unit Commitment problem using a Mixed Integer Linear Programming approach in MATLAB. I was able to do this with some help from this link: https://www.mathworks.com/examples/matlab/community/36286-script-9-unit-commitment. I just had to make some modifications since I wanted a piecewise linear formulation of the cost curve. This is the modified code:

%Get the Load Demand
Load = xlsread('Load_Demand.xlsx', -1);

%Get the Number of Hours to perform optimization operations
nHours = numel(Load);
Time = (1:nHours);

%Enter the Number of Generating Units in the System
num_of_gen = str2double(inputdlg('Enter the Number of Units in the System','Number of Generating Units'));


%Enter Generator Data from Excel File
%Create Matrix to store data
gen_data = zeros(num_of_gen, 8);

gen_data = xlsread('Gen_Data.xlsx', -1);

%Pull apart generator properties

%Get Start Up Cost data
StartupCost  = (gen_data(:, 8))';

%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, 1) * gen_data(num, 5) * gen_data(num, 5)) + (gen_data(num, 2) * gen_data(num, 5)) + gen_data(num, 3);

end

%Get the Minimum up and down times data
MinimumUpTime = (gen_data(:, 6))';
MinimumDownTime = (gen_data(:, 7))';

%==========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, 4) - gen_data(num, 5)) / 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

%Get Maximum generation level for the actual units
MaxPowerGenLevel = gen_data(:, 4)';

%Set Minimum generation level to 0
MinGenerationLevel = zeros(1, (num_of_gen*K));

%Get Minimum 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, 5);
            
        elseif k == (K + 1)
            
            Power_Segment(num, k) = gen_data(num, 4);
            
        else
            
            Power_Segment(num, k) = gen_data(num, 5) + ((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, 1) * Power_Segment(num, k) * Power_Segment(num, k)) + (gen_data(num, 2) * Power_Segment(num, k)) + gen_data(num, 3);
        
    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;

idxHr1 = 1;
idxHr2ToEnd=2:nHours;

maxGenConst = repmat(MaxGenerationLevel,nHours,1);
minGenConst = repmat(MinGenerationLevel,nHours,1);
minPowConst = repmat(MinPowerGenLevel,nHours,1);
maxPowConst = repmat(MaxPowerGenLevel,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(nHours, num_of_gen);

for num = 1:1:num_of_gen
    
    sub_maxGenConst(:, num) = maxGenConst(:, (num * K));
    
end


new_maxGenConst = sub_maxGenConst.*isOn;


%sub_new_maxGenConst = zeros(24, (num_of_gen*K));
%Initialise sub_new_maxGenConst
sub_new_maxGenConst = repmat(new_maxGenConst(:, 1), 1, K);
for num = 2:1:num_of_gen
     
    sub_new_maxGenConst = [sub_new_maxGenConst repmat(new_maxGenConst(:, num), 1, K)];
    
end


%minGenConst 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 minGenConst to a 24x(num_of_gen) matrix
%And then transforms it back to a 24x(num_of_gen*K) optimization matrix
sub_minGenConst = zeros(nHours, num_of_gen);

for num = 1:1:num_of_gen
    
    sub_minGenConst(:, num) = minGenConst(:, (num * K));
    
end


new_minGenConst = sub_minGenConst.*isOn;


%sub_new_minGenConst = zeros(24, (num_of_gen*K));
%Initialise sub_new_minGenConst
sub_new_minGenConst = repmat(new_minGenConst(:, 1), 1, K);
for num = 2:1:num_of_gen
     
    sub_new_minGenConst = [sub_new_minGenConst repmat(new_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);


%Spinning Reserve Constraint
powerprob.Constraints.spinningReserve = sum((maxPowConst.*isOn),2) >= Load * 1.10;


%Upper Bound Constraints which change with isOn
%Only gen power when plant is on
%If isOn=0 power must be zero, if isOn=1 power must be less than maxGenConst
powerprob.Constraints.powerOnlyWhenOn = power <= sub_new_maxGenConst; 

%Lower Bound Constraints which change with isOn
% if on, meet MinGenerationLevel
% if isOn=0 power >= 0, if isOn=1 power must be more than minGenConst
powerprob.Constraints.meetMinGenLevel = power >= sub_new_minGenConst; 

%enforce startup=1 when moving from off to on for the first hour
powerprob.Constraints.startupHr1 = isOn(idxHr1,:) - startup(idxHr1,:) <= 0;

%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(nHours, 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:nHours
    
    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(nHours, num_of_gen);

for h = 1:1:nHours
    
    for num = 1:1:num_of_gen
        
        sub_Hourly_Cost(h, num) = ((gen_data(num, 1) * Total_Power(h, num) * Total_Power(h, num)) + (gen_data(num, 2) * Total_Power(h, num)) + gen_data(num, 3));
        
    end
    
end

sub_Hourly_Cost = sub_Hourly_Cost .* sol.isOn;
        
%Find Total Cost of Operation
Hourly_Cost = zeros(nHours, 1);

for h = 1:1:nHours

    for num = 1:1:num_of_gen

        Hourly_Cost(h, 1) = Hourly_Cost(h, 1) + sub_Hourly_Cost(h, num);

    end

end

Now I need to include some energy storage devices in the system which can act as a load when charging and a generator when discharging. However, the manner in which the energy storage is used is dependent on whether a new generating unit is required when moving from one hour to the next. Thus, my problem lies with determining whether a new unit is required before actually determining how to use the energy storage, and then when I determine this I need to modify the load demand (if charging, load will increase by charging rate and if dicharging load will decrease by discharging rate) and re-run the simulation to meet the new load demand.

One way I am thinking of doing this is by removing the following constraint:

powerprob.Constraints.isDemandMet = sum(power,2) >= Load - sum((minPowConst.*isOn),2);

and include the following. This way I hope to run the simulation without considering the energy storage and after determining how to use the energy storage, re-run the simulation. If this does not make sense, I would appreciate some guidance in this matter, as to how I can properly implement my idea.

for kk = 1:nHours
    if kk == 1
        %Energy Storage is Not operational for first hour
        BESS_is_Load = 0;
        BESS_is_Gen = 0;
    else
        jj = 1:num_of_gen;
        powerprob.Constraints.isDemandMet1 = sum(power,2) >= Load - sum((minPowConst.*isOn),2);
        %Determine if New Unit is Required
        if sum(isOn(kk, jj)) > sum(isOn((kk-1), jj))
            %Determine if Storage is Empty

            %Update BESS_is_Load and BESS_is_Gen values
            
        else
            %Determine if Online Units are operating at MAX capacity

            %Update BESS_is_Load and BESS_is_Gen values

        end 
    end
    powerprob.Constraints.isDemandMet2 = sum(power,2) >= Load - sum((minPowConst.*isOn),2) + (BESS_is_Load * Charging_Rate) - (BESS_is_Gen * Discharging_Rate);
end       

[ayayda amira]

do you have the unit Commitement Data for this example.I will be very grateful to you it if you can send me the file to: [e-mail].

[Rehman Zafar]

Hi, it is very interesting could you please send me the data file that you have used in this code at [e-mail] Thank you!

[Priya Nayar]

Hello Micah,

Is battery size and cost minimisation possible with an extension tothis code, if so please share your support on [e-mail]

Thanks

[Niccolò Callarà]

Hello Micah,

this is just the case I'm trying to analyze and study these days, I think this code will help me a lot. I will be very grateful if you send this code to my email ([e-mail]).

Thank you so much!

[BABA SA]

dear Sir,

i am studying unit commitment as at now. Can you please help me with the codes please. It will be a great help.

Email:[e-mail]

thanks

[Zeeshan Hameed]

Sir Can you please send the data file for this program. That will be very much helpful for me. Thank you sir My email: [e-mail]

[Hussain Kakar]

sir can you please send me this file my email is [e-mail]

[Jennifer Solorzano]

Sir can you please send me this file my email is [e-mail]. Thank you Sir