Solving the Unit Commitment Problem using optimproblem

I am trying to solve the Unit Commitment problem using the optimproblem function. I will explain my problem with an example. Lets say I have three generators which will be operated to meet various loads for a 24 hours period. Now the cost curves for these units are quadratic, so piecewise linearization is required if I am going to use the intlinprog function. I have successfully linearized the system and the problem now involves minimizing:

     Total C(P) = C_1(P_1min)T_1 + s_11*P_11*T_1 + s_12*P_12*T_1 + … + s_1n*P_1n*T_1
            + C_2(P_2min)T_2 + s_21*P_21*T_2 + s_22*P_22*T_2 + … + s_2n*P_2n*T_2
            + C_3(P_3min)T_3 + s_31*P_31*T_3 + s_32*P_32*T_3 + … + s_3n*P_3n*T_3

Where Sik is the slope for each segment of the linearized curve, Pik are the power increments (which I need to determine) for each segment. T_i is a binary variable which keeps track of the state of unit and Ci(Pimin) is the cost at Pimin. Now based on my research, the problem should now be formulated as minimize:

    Total C(P) = C_1(P_1min)T_1 + s_11*Z_11 + s_12*Z_12 + … + s_1n*Z_1n
            + C_2(P_2min)T_2 + s_21*Z_21 + s_22*Z_22 + … + s_2n*Z_2n
            + C_3(P_3min)T_3 + s_31*Z_31 + s_32*Z_32 + … + s_3n*Z_3n

Where Ti = 0 implies that Zik = 0 and Ti = 1 implies that Zik = Pik.

Now I have two optimization variable:

1. Pik
2. Ti

With the optimization problem being made up of two part:

1. Sum of ALL Sik*Pik values
2. Sum of Ci(Pimin)*Ti values.

The following code illustrates this (with the addition of another optimization variable which isn't so important):

 %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);

    %Costs
powerCost = sum(power*f,1);
isOnCost = sum(isOn*OperatingCost',1);
startupCost = sum(startup*StartupCost',1);

 % set objective
powerprob.Objective = powerCost + isOnCost + startupCost;

   % power generation over all plans meets hourly demand
powerprob.Constraints.isDemandMet = sum(power,2) >= Load;

    % 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 <= maxGenConst.*isOn;   

    % if on, meet MinGenerationLevel
% if isOn=0 power >= 0, if isOn=1 power must be more than minGenConst
powerprob.Constraints.meetMinGenLevel = power >= maxGenConst.*isOn;

My problem lies in this calculation: maxGenConst.*isOn. These have different sizes as maxGenConst is a nHours * (num_of_gen*num_of_segments) whereas isOn is a nHour * num_of_gen. The reason this occurs is because after linearization, the cost equation becomes (num_of_segments + 1) equations. I would appreciate some help in formulating the problem so I can use the optimproblem function. I have attached my full code and data.

[Matt J]

My problem lies in this calculation: maxGenConst.*isOn. These have different sizes as maxGenConst is a nHours * (num_of_gen*num_of_segments) whereas isOn is a nHour * num_of_gen

If they have different sizes, then in what way does it make sense to try multiplying them together?

[Micah Mungal]

I attempted to use another approach to solve the unit commitment problem using the optimproblem function. Did not know this function existed until a couple days ago. But I like that I can just state the constraints using this approach rather than formulating it myself as that way is more prone to humam error.

This multiplication is suppose to be a piecewise multiplication to update the upper limit of maxgenerationlimit based on whether the unit is on or off. But based on the way I formulated my objective function and due to the many segments of the linearized curve, I ended up with the variables being of different size. (A group of segments actually make up the cost function for one generator). So if I have three generator and I used 4 segments to linearize the cost for each generator I end up with 3 x 4 = 12 power segments in total. Whereas the isOn should only be used to keep track of the 3 generators and not all 12 segments.

I attempted to index values in the isOn variable to make it the same size as the maxGenConst variable but ran into alot of problems. I dont think you can access the optimization variable. I would not have run into the problem if I used the quadratic formulation for the cost curve but MATLAB doed not have a MIQP solver so I am stuck with linearizing the cost function. I am just looking for some advice as to what I should do.

[Micah Mungal]

I was able to solve my problem. I made maxGenConst the same size as isOn. Then I did the dot multiplication with isOn and was able to index the optimization expression and reshape into the proper dimensions. I have attached my code with data for a test system for anyone who will need it. This code considers power generation constraints, minimum up and down times, start up costs and the basic load demand constraints. I have included the data for a 6 unit system with the necessary data and the load demand over a 24hr period. This code does not consider losses in the system.

[Sivateja Maturu]

Can you please explain how you have defined the values of a, b, c for the power plants defined in your Gen_Data?

Answer 1

[Micah Mungal]

I was able to solve my problem. I made maxGenConst the same size as isOn. Then I did the dot multiplication with isOn and was able to index the optimization expression and reshape into the proper dimensions. I have attached my code with data for a test system for anyone who will need it. This code considers power generation constraints, minimum up and down times, start up costs and the basic load demand constraints. I have included the data for a 6 unit system with the necessary data and the load demand over a 24hr period. This code does not consider losses in the system.

[MUTHUSELVI G]

In this code Error in UC (line 178) powerprob = optimvar; is coming . what will i do?

[Micah Mungal]

What do you mean by "is coming"?

[ali gh]

An error appears after running the code

Subscripted assignment dimension mismatch.

Error in Untitled2 (line 244) sub_maxGenConst(:, num) = maxGenConst(:, (num * K));

[Martin Cruz]

tiene errores tu huevada

[Kadyrbek]

Hello. Could you please send the code. [e-mail]

[Walter Roberson]

@MUTHUSELVI G

The code needs Optimization Toolbox, R2017b or later.

Answer 2

[amrit singh]

Undefined function or variable 'optimproblem'.

Error in ELD_LP (line 168)

powerprob = optimproblem;

[Walter Roberson]

The code requires Optimization Toolbox, R2017b or later.

Answer 3

[praveen kumar]

%==========================================================================
% ECONOMIC LOAD DISPATCH - DIAGNOSTIC AND WORKING VERSION
%==========================================================================
clc;
clear;
close all;
%% ==================== DATA INPUT ====================
Load = [140; 166; 180; 196; 220; 240; 267; 283.4; 308; 323; 340; 350; ...
    300; 267; 220; 196; 220; 240; 267; 300; 267; 235; 196; 166];
nHours = numel(Load);
% Generator Data: [b, a, c, Pmax, Pmin, MinUpTime, MinDownTime, StartupCost]
gen_data = [2       200     0.00375    200    50     1    1    176;
    1.75    257     0.0175     80     20     2    2    187;
    1       300     0.0625     40     15     1    1    113;
    3.25    400     0.00834    35     10     1    2    267;
    3       515     0.025      30     10     2    1    180;
    3       515     0.05       25     12     1    1    113];
num_of_gen = size(gen_data, 1);
b_coeff = gen_data(:, 1);
a_coeff = gen_data(:, 2);
c_coeff = gen_data(:, 3);
Pmax = gen_data(:, 4);
Pmin = gen_data(:, 5);
fprintf('========== FEASIBILITY CHECK ==========\n');

========== FEASIBILITY CHECK ==========

fprintf('Peak Load: %.2f MW\n', max(Load));

Peak Load: 350.00 MW

fprintf('Total Capacity: %.2f MW\n', sum(Pmax));

Total Capacity: 410.00 MW

fprintf('Min Load: %.2f MW\n', min(Load));

Min Load: 140.00 MW

fprintf('Total Pmin (all ON): %.2f MW\n\n', sum(Pmin));

Total Pmin (all ON): 117.00 MW

% THE ISSUE: Min load (140) is greater than what some generators can provide alone
% but less than total Pmin if all units are ON (117 MW is fine)
% The problem is the PIECEWISE LINEARIZATION creating infeasibility
%% ========== SIMPLIFIED APPROACH WITHOUT PIECEWISE ==========
% Use quadratic costs directly - simpler and more reliable
fprintf('========== USING SIMPLIFIED FORMULATION ==========\n');

========== USING SIMPLIFIED FORMULATION ==========

% Create optimization problem
prob = optimproblem('ObjectiveSense', 'minimize');
% Decision Variables
% P(t,g) = power output of generator g at hour t
P = optimvar('P', nHours, num_of_gen, 'LowerBound', 0);
% u(t,g) = 1 if generator g is ON at hour t
u = optimvar('u', nHours, num_of_gen, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
% v(t,g) = 1 if generator g starts up at hour t
v = optimvar('v', nHours, num_of_gen, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
fprintf('Variables: %d continuous, %d binary\n', nHours*num_of_gen, 2*nHours*num_of_gen);

Variables: 144 continuous, 288 binary

%% ========== OBJECTIVE FUNCTION ==========
% For optimization, we'll use linear approximation of costs
% Cost ≈ marginal cost * power + fixed cost when on
% Calculate marginal cost at mid-point
marginal_cost = zeros(num_of_gen, 1);
fixed_cost = zeros(num_of_gen, 1);
startup_cost = gen_data(:, 8);
for g = 1:num_of_gen
    P_mid = (Pmin(g) + Pmax(g)) / 2;
    marginal_cost(g) = 2 * a_coeff(g) * P_mid + b_coeff(g);
    fixed_cost(g) = a_coeff(g) * Pmin(g)^2 + b_coeff(g) * Pmin(g) + c_coeff(g);
end
% Objective: minimize total cost
fuel_cost = sum(P * marginal_cost, 'all');
operating_cost = sum(u * fixed_cost, 'all');
startup_cost_total = sum(v * startup_cost, 'all');
prob.Objective = fuel_cost + operating_cost + startup_cost_total;
fprintf('Objective function defined\n');

Objective function defined

%% ========== CONSTRAINTS ==========
fprintf('Adding constraints...\n');

Adding constraints...

% 1. Meet demand
prob.Constraints.demand = sum(P, 2) >= Load;
fprintf('  Demand: %d constraints\n', nHours);

  Demand: 24 constraints

% 2. Minimum generation when ON
for g = 1:num_of_gen
    prob.Constraints.(['min_' num2str(g)]) = P(:, g) >= Pmin(g) * u(:, g);
end
fprintf('  Min generation: %d constraints\n', num_of_gen * nHours);

  Min generation: 144 constraints

% 3. Maximum generation when ON  
for g = 1:num_of_gen
    prob.Constraints.(['max_' num2str(g)]) = P(:, g) <= Pmax(g) * u(:, g);
end
fprintf('  Max generation: %d constraints\n', num_of_gen * nHours);

  Max generation: 144 constraints

% 4. Startup definition (relaxed)
prob.Constraints.startup = v(2:end, :) >= u(2:end, :) - u(1:end-1, :);
fprintf('  Startup: %d constraints\n', (nHours-1) * num_of_gen);

  Startup: 138 constraints

fprintf('Total constraints: %d\n', nHours + 2*num_of_gen*nHours + (nHours-1)*num_of_gen);

Total constraints: 450

%% ========== SOLVE ==========
fprintf('\n========== SOLVING ==========\n');

========== SOLVING ==========

options = optimoptions('intlinprog', ...
    'MaxTime', 600, ...
    'Display', 'iter', ...
    'RelativeGapTolerance', 0.1, ...
    'MaxNodes', 1e6);
fprintf('Starting solver...\n\n');

Starting solver...

tic;
[sol, fval, exitflag, output] = solve(prob, 'Options', options);

Solving problem using intlinprog.
Running HiGHS 1.7.1: Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 2e+02]
  Cost   [1e+02, 5e+05]
  Bound  [1e+00, 1e+00]
  RHS    [1e+02, 4e+02]
Presolving model
450 rows, 426 cols, 1134 nonzeros  0s
360 rows, 364 cols, 972 nonzeros  0s
299 rows, 333 cols, 851 nonzeros  0s

Solving MIP model with:
   299 rows
   333 cols (244 binary, 0 integer, 10 implied int., 79 continuous)
   851 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   116666184.2554  inf                  inf        0      0      0         0     0.0s
         0       0         0   0.00%   175071916.2681  inf                  inf        0      0      2       100     0.0s
 R       0       0         0   0.00%   175620875.9558  175706768.0089     0.05%       52     13      2       117     0.0s

Solving report
  Status            Optimal
  Primal bound      175706768.009
  Dual bound        175638839.359
  Gap               0.0387% (tolerance: 10%)
  Solution status   feasible
                    175706768.009 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.01 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             1
  LP iterations     118 (total)
                    0 (strong br.)
                    17 (separation)
                    0 (heuristics)

Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.RelativeGapTolerance = 0.1. The intcon variables are integer
within tolerance, options.ConstraintTolerance = 1e-06.

solve_time = toc;
%% ========== DISPLAY RESULTS ==========
fprintf('\n========== RESULTS ==========\n');

========== RESULTS ==========

fprintf('Exit flag: %d\n', exitflag);

Exit flag: 1

fprintf('Solve time: %.2f seconds\n', solve_time);

Solve time: 0.89 seconds

if exitflag >= 0
    fprintf('STATUS: SUCCESS!\n');
    fprintf('Total Cost: $%.2f\n', fval);
    
    P_sol = sol.P;
    u_sol = sol.u;
    v_sol = sol.v;
    
    % Round binary variables
    u_sol = round(u_sol);
    v_sol = round(v_sol);
    
    fprintf('\n========== GENERATION SCHEDULE ==========\n');
    fprintf('Hour\tLoad\tTotal\tUnits\t  G1\t  G2\t  G3\t  G4\t  G5\t  G6\n');
    fprintf('----\t----\t-----\t-----\t----\t----\t----\t----\t----\t----\n');
    
    for t = 1:nHours
        fprintf('%2d\t%4.0f\t%5.1f\t  %d\t', t, Load(t), sum(P_sol(t,:)), sum(u_sol(t,:)));
        for g = 1:num_of_gen
            if u_sol(t,g) > 0.5
                fprintf('%5.1f\t', P_sol(t,g));
            else
                fprintf('  -  \t');
            end
        end
        fprintf('\n');
    end
    
    % Calculate actual costs
    actual_gen_cost = 0;
    actual_op_cost = 0;
    for t = 1:nHours
        for g = 1:num_of_gen
            if u_sol(t,g) > 0.5
                P_val = P_sol(t,g);
                actual_gen_cost = actual_gen_cost + a_coeff(g)*P_val^2 + b_coeff(g)*P_val + c_coeff(g);
            end
        end
    end
    actual_startup_cost = sum(sum(v_sol)) * mean(startup_cost);
    
    fprintf('\n========== COST BREAKDOWN ==========\n');
    fprintf('Generation Cost: $%.2f\n', actual_gen_cost);
    fprintf('Startup Cost:    $%.2f (approx)\n', actual_startup_cost);
    fprintf('Total:           $%.2f (approx)\n', actual_gen_cost + actual_startup_cost);
    
    fprintf('\n========== STARTUPS ==========\n');
    for g = 1:num_of_gen
        fprintf('Generator %d: %d startups\n', g, sum(v_sol(:,g)));
    end
    
    % Verify demand is met
    fprintf('\n========== VERIFICATION ==========\n');
    for t = 1:nHours
        total_gen = sum(P_sol(t,:));
        if total_gen < Load(t) - 0.01
            fprintf('WARNING: Hour %d - Generation %.2f < Load %.2f\n', t, total_gen, Load(t));
        end
    end
    fprintf('All demands satisfied!\n');
    
    % Plot results
    figure('Position', [100, 100, 1400, 700]);
    
    % Subplot 1: Stacked generation
    subplot(2,2,1);
    area(1:nHours, P_sol);
    hold on;
    plot(1:nHours, Load, 'k-', 'LineWidth', 2);
    xlabel('Hour');
    ylabel('Power (MW)');
    title('Generation Dispatch');
    legend('G1','G2','G3','G4','G5','G6','Load','Location','eastoutside');
    grid on;
    
    % Subplot 2: Unit commitment
    subplot(2,2,2);
    imagesc(u_sol');
    colormap([1 1 1; 0 0.5 0]);
    xlabel('Hour');
    ylabel('Generator');
    title('Unit Commitment Status');
    yticks(1:num_of_gen);
    yticklabels({'G1','G2','G3','G4','G5','G6'});
    colorbar('Ticks',[0.25 0.75],'TickLabels',{'OFF','ON'});
    
    % Subplot 3: Load vs Generation
    subplot(2,2,3);
    plot(1:nHours, Load, 'b-o', 'LineWidth', 2, 'MarkerFaceColor', 'b');
    hold on;
    plot(1:nHours, sum(P_sol,2), 'r-s', 'LineWidth', 2, 'MarkerFaceColor', 'r');
    xlabel('Hour');
    ylabel('Power (MW)');
    title('Load vs Total Generation');
    legend('Load','Generation','Location','best');
    grid on;
    
    % Subplot 4: Generator utilization
    subplot(2,2,4);
    utilization = zeros(num_of_gen, 1);
    for g = 1:num_of_gen
        hours_on = sum(u_sol(:,g));
        avg_output = sum(P_sol(:,g)) / hours_on;
        utilization(g) = avg_output / Pmax(g) * 100;
    end
    bar(utilization);
    xlabel('Generator');
    ylabel('Average Utilization (%)');
    title('Generator Utilization When ON');
    xticklabels({'G1','G2','G3','G4','G5','G6'});
    grid on;
    
    sgtitle('Economic Load Dispatch Results', 'FontWeight', 'bold', 'FontSize', 14);
    
else
    fprintf('STATUS: FAILED\n');
    fprintf('Message: %s\n', output.message);
    
    % Additional diagnostics
    fprintf('\n========== DIAGNOSTICS ==========\n');
    fprintf('The problem is infeasible. Possible reasons:\n');
    fprintf('1. Check if any hour load > total capacity\n');
    for t = 1:nHours
        if Load(t) > sum(Pmax)
            fprintf('   Hour %d: Load %.0f > Capacity %.0f\n', t, Load(t), sum(Pmax));
        end
    end
    fprintf('2. All capacity checks passed - issue is with constraints\n');
end

STATUS: SUCCESS!

Total Cost: $175706768.01

========== GENERATION SCHEDULE ==========

Hour	Load	Total	Units	  G1	  G2	  G3	  G4	  G5	  G6

----	----	-----	-----	----	----	----	----	----	----

 1	 140	140.0	  4	

  -  	

 40.0	 40.0	 35.0	

  -  	

 25.0	

 2	 166	166.0	  5	

  -  	

 36.0	 40.0	 35.0	 30.0	 25.0	

 3	 180	180.0	  5	

  -  	

 50.0	 40.0	 35.0	 30.0	 25.0	

 4	 196	196.0	  5	

  -  	

 66.0	 40.0	 35.0	 30.0	 25.0	

 5	 220	220.0	  6	

 50.0	 40.0	 40.0	 35.0	 30.0	 25.0	

 6	 240	240.0	  6	

 50.0	 60.0	 40.0	 35.0	 30.0	 25.0	

 7	 267	267.0	  6	

 57.0	 80.0	 40.0	 35.0	 30.0	 25.0	

 8	 283	283.4	  6	

 73.4	 80.0	 40.0	 35.0	 30.0	 25.0	

 9	 308	308.0	  6	

 98.0	 80.0	 40.0	 35.0	 30.0	 25.0	

10	 323	323.0	  6	

113.0	 80.0	 40.0	 35.0	 30.0	 25.0	

11	 340	340.0	  6	

130.0	 80.0	 40.0	 35.0	 30.0	 25.0	

12	 350	350.0	  6	

140.0	 80.0	 40.0	 35.0	 30.0	 25.0	

13	 300	300.0	  6	

 90.0	 80.0	 40.0	 35.0	 30.0	 25.0	

14	 267	267.0	  6	

 57.0	 80.0	 40.0	 35.0	 30.0	 25.0	

15	 220	220.0	  6	

 50.0	 40.0	 40.0	 35.0	 30.0	 25.0	

16	 196	196.0	  5	

  -  	

 66.0	 40.0	 35.0	 30.0	 25.0	

17	 220	220.0	  6	

 50.0	 40.0	 40.0	 35.0	 30.0	 25.0	

18	 240	240.0	  6	

 50.0	 60.0	 40.0	 35.0	 30.0	 25.0	

19	 267	267.0	  6	

 57.0	 80.0	 40.0	 35.0	 30.0	 25.0	

20	 300	300.0	  6	

 90.0	 80.0	 40.0	 35.0	 30.0	 25.0	

21	 267	267.0	  6	

 57.0	 80.0	 40.0	 35.0	 30.0	 25.0	

22	 235	235.0	  6	

 50.0	 55.0	 40.0	 35.0	 30.0	 25.0	

23	 196	196.0	  5	

  -  	

 66.0	 40.0	 35.0	 30.0	 25.0	

24	 166	166.0	  5	

  -  	

 36.0	 40.0	 35.0	 30.0	 25.0	

========== COST BREAKDOWN ==========

Generation Cost: $90402132.01

Startup Cost:    $690.67 (approx)

Total:           $90402822.68 (approx)

========== STARTUPS ==========

Generator 1: 2 startups
Generator 2: 1 startups
Generator 3: 0 startups
Generator 4: 0 startups
Generator 5: 1 startups
Generator 6: 0 startups

========== VERIFICATION ==========

All demands satisfied!

fprintf('\nDone.\n');

Done.