Solving the Unit Commitment Problem using optimproblem

Question

Micah Mungal 2018-5-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/403089-solving-the-unit-commitment-problem-using-optimproblem

评论： Walter Roberson 2025-10-28，2:28

UC Code.zip

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:

Pik
Ti

With the optimization problem being made up of two part:

Sum of ALL Sik*Pik values
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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Micah Mungal 2018-5-30

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 2018-6-11

UC using MILP.zip

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 2021-6-16

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 2018-6-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/403089-solving-the-unit-commitment-problem-using-optimproblem#answer_324152

UC using MILP.zip

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Kadyrbek 2022-9-26

Hello. Could you please send the code. k.kozhobekov＠yahoo.com

Walter Roberson 2025-10-28，2:28

@MUTHUSELVI G

The code needs Optimization Toolbox, R2017b or later.

请先登录，再进行评论。

Answer 2

amrit singh 2019-5-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/403089-solving-the-unit-commitment-problem-using-optimproblem#answer_375718

Undefined function or variable 'optimproblem'.

Error in ELD_LP (line 168)

powerprob = optimproblem;

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2025-10-28，1:46

The code requires Optimization Toolbox, R2017b or later.

请先登录，再进行评论。

Answer 3

praveen kumar 2025-10-28，1:34

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/403089-solving-the-unit-commitment-problem-using-optimproblem#answer_1571454

编辑：Walter Roberson 2025-10-28，1:47

在 MATLAB Online 中打开

%==========================================================================

% 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

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

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

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

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Solving the Unit Commitment Problem using optimproblem

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（2 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Solving the Unit Commitment Problem using optimproblem

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（2 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论