Using intlinprog for Unit Commitment Problem

I am trying to solve the Unit Commitment problem using intlinprog. Now I have already solved the Economic Load Dispatch problem using linprog and I intend on expanding the problem to solve the Unit Commitment problem.

For the economic load dispatch problem, I have three generating units which are assumed to be committed and I need to determine how much each power each unit must produce to meet the load demand at the lowest cost. Assume that the cost function for a unit is given as:

C(P) = a*P^2 + b*P + c.

I was able to linearize this cost function in the form of:

C(P) = K(P_1, P_2,...,P_n) = C(P_min) + s_1*P_1 + s_2*P_2 + … + s_n*P_n

where s_1, s_2, … , s_n are the slopes of the individual segments.

Now P is originally bounded by an upper and lower bound, and after linearization P_1, P_2, …. ,P_n are now bounded by values based on the originally bounds. Also the sum of all P's must meet a load demand. I was able to successfully formulate this problem in MATLAB and solve it using the linprog function. See attached files.

However, I now have to expand this concept to the Unit Commitment problem, where I have to now determine which units should be committed in the first place and then solve for the economic load dispatch problem. So the problem now involves a variable u = 0, or u = 1 which determines whether the unit is on/off. The objective function now becomes minimize (for three units):

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

And based on research, I determined the problem must be formulated as:

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

where u_i = 0, implies that Z_ik = 0 and u_i = 1, implies that Z_ik = Pik. So from this it is simple to formulate the objective function and the unknow variables will now be the u_i and Z_ik variables.

My problem lies in creating the upper bounds. The lower bound will always be zero but the upper bound will change depending on the value of u_i. How do I formulate this problem properly?

Answer 1

[Matt J]

I think you should formulate it this way

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

with the constraints

   0<= P_ij <=Upper_ij

   T_i >= sum_j(P_ij)/ sum_j(Upper_ij)

   T_i binary

in addition to whatever minimum output sum(P_ij) must satisfy.

This works because if T_i=0 then all corresponding P_ij must equal zero by virtue of the first two constraints above. Conversely if even one P_ij>0 then likewise T_i>0 which in turn forces T_i=1 because of the integer/binary constraint on T_i.

[Micah Mungal]

Thanks for the response. Just making sure I understanding what everything. If I formulate the problem the way you stated then my x (with respect to the way intlinprog expects the problem to formulated) should be:

x = [T_1; P_11; P_12;..., T_2; P21; P22;....T_n; P_n1, ....P_nn);

For the load constraint where the sum of all power generated should equal to the load, then Aeq should following something like this:

Aeq = [P_1min 1 1.....P_2min 1 1 .....P_nmin 1 1....] and beq = load

Now this is the part I am unsure about. You said to include a constraint:

T_i >= sum_j(P_ij)/ sum_j(Upper_ij)

Therefore to use this with the intinprog, I determined that for the inequality constraint Ax >= b

it must be rewritten as follows:

sum_j(Upper_ij)*T_i - sum_j(P_ij) >= 0

Therefore

A = [sum_j -1 -1 -1 .....sum_j -1 -1 -1...] with b = [0 0 0 .... 0]

This formulation does not seem right too me.

[Matt J]

The only thing wrong that I see is that it should be A*x<=b, so the A should be

    A = [-sum_j +1 +1 +1 ...]

[Micah Mungal]

Consider this T_1 >= sum_j(P_1j)/ sum_j(Upper_1j)

T_2 >= sum_j(P_2j)/ sum_j(Upper_2j)

T_3 >= sum_j(P_3j)/ sum_j(Upper_3j)

Now as it is, x is a (ij + i) x 1 matrix, since for each generator there will be j terms plus the C_i(P_imin) value for each generator. If I use this formulation for x, then I cant use A = [sum_j -1 -1 -1 .....sum_j -1 -1 -1...]

as this does not consider the contraint T_i >= sum_j(P_ij)/ sum_j(Upper_ij) separately but rather as one big sum.

[Matt J]

That is correct. A must have N rows, one for each generator, and each row must contain non-zero coefficients only where applicable to that generator:

   A = [-sum_Upper_1j +1 +1 +1                0  0  0  0  0 0 0 0 .... ; 
                    0  0  0  0    -sum_Upper_2j +1 +1 +1  0 0 0 0 ...  ;
         ...]

[Micah Mungal]

Hey man, thanks in a million. Adding the constraint, T_i >= sum_j(P_ij)/ sum_j(Upper_ij) worked beautifully.

I have one addition constraint to add. Not sure if you can help me out with this one. Now that I have built the code, it performs the simulation for 1 time interval. Now I can adjust the code to include different load values for 24 time interval (to represent a day). However, the result would not be realistic as there is a minimum up and down time constraint associated with each unit. For example if a unit is turned on (i.e. T_i = 1) at a specific time interval, it must remain on for a certain number of hours and in this case a number of time interval before it can be turned off again (T_i if 1, must be 1 for some time before it can be T_i = 0). Similarly, if a unit is turn off (i.e. T_i = 0) at a specific interval it must remain off for a specific time interval before it can be turned on again.

Now I am not sure how I should go about integrating this with the intlinprog function. Should I do some kind of if statement and do a check at each interval or should it be a constraint that is included when the intlinprog function is called?

[Alan Weiss]

Take a look at this example. I think that the ideas it shows about how to charge for turning on and off a generator will generalize to your requirement.

Alan Weiss

MATLAB mathematical toolbox documentation

[Matt J]

Hey man, thanks in a million. Adding the constraint, T_i >= sum_j(P_ij)/ sum_j(Upper_ij) worked beautifully.

You're welcome, but please Accept-click the answer to certify that your original question was addressed.

I have one addition constraint to add.

If Alan's suggestion doesn't help, I would post this as a new question.

[Micah Mungal]

I had a good look at Alan's suggestion but it doesn't help.

[Micah Mungal]

My approach to this problem currently involves using the intlinprog to solve the problem one hour at a time while keeping track of the minimum up and down time and using if statements to set the upper bounds for binary variable (used to specify if a units is on/off) to 1 or 0 depending on the outcome of the if statement. For example, lets say the unit if turned on in the 3rd hour and has a minimum up time of 2 hours, then I will set both the upper and lower bound for the binary variable as 1 thereby restricting the unit from turning off.

[Matt J]

After some thought, I think you can do it by including inequality constraints on the T_i differences. Suppose for example that the minimum up time is 3 hours. Then you would impose for all j,

T_j >= T_{j-1}-T_{j-2}  % forces off-to-on transition to hold for 3 hours
T_j >= T_{j-2}-T_{j-3}
T_j >= T_{j-3}-T_{j-4}

And, if the minimum off time is 4 hours for example, you impose, for all k,

T_k <= T_{k-1}-T_{k-2}  + 1; % forces on-to-off transition to hold for 4 hours
T_k <= T_{k-2}-T_{k-3}  + 1;
T_k <= T_{k-3}-T_{k-4}  + 1;
T_k <= T_{k-4}-T_{k-5}  + 1;

[Micah Mungal]

Ok, after considering this approach and assuming that j and k in the above formulation is the current hour, I have concluded that I will need to reformulate the unknown x matrix to include all 24 hours. Is this correct. I can formulate x to take the form of for three units:

[T_11  T_21  T_31  … T_(24)1
 P_111 P_211 P_311 … P_(24)11
 P_112 P_212 P_312 … P_(24)12
    .
    .
 T_12  T_22  T_32  … T_(24)2
 P_121 P_221 P_321 … P_(24)21
 P_122 P_222 P_322 … P_(24)22
    .
    .
 T_13  T_23  T_33  … T_(24)3
 P_131 P_231 P_321 … P_(24)31
 P_132 P_232 P_322 … P_(24)32
   .
 P_13n P_23n P_33n … P_(24)3n
]

Where the first subscript indicates the hour, the second indicates the unit number and wrt to the power variable P, the last subscript indicates the segment after linearization.

That is x becomes are (number of units*(K+1)) X 24 matrix where K is the number of segments used to linearize the curve and 24 for the 24 hours considered.

Now A_eq, b_eq, the lb and ub are simple to formulate. Also for the inequality constraint A and b based on the following constraint is also simple to formulate.

   0<= P_ij <=Upper_ij
   T_i >= sum_j(P_ij)/ sum_j(Upper_ij)
   T_i binary

However, I am having some difficulty including the minimum up and down time constraint. Since intlinprog only allows for one inequality constraint of A.x <= b, then I must find some way of including it in the above constraint. Well this is what I am thinking. I could be wrong.

[Matt J]

Also for the inequality constraint A and b based on the following constraint is also simple to formulate.

Note that the constraint 0<= P_ij <=Upper_ij should be handled with lb,ub, not with A,b.

Since intlinprog only allows for one inequality constraint of A.x <= b, then I must find some way of including it in the above constraint.

I'm not sure what you mean by "only allows for one". A and b are matrices with N rows, each of which defines a constraint. Just include a row for every inequality constraint that you need.

[Micah Mungal]

So I tried to formulate the problem as stated in my last post on the 5th June and realised that I cannot formulate x as stated. x must remain as a (num of gen * (K + 1)) X 1 matrix. The K + 1 incorporates the binary variable which indicates whether the unit is ON/OFF. It is f (based on the documentation for intlinprog) that now becomes a (24 * num of gen * (K + 1)) X 1 matrix. The 24 indicates 24 hours. The number of unknown variable now becomes (24 * num of gen * (K + 1)).

I then formulated Aeq as a 24 X (24 * num of gen * (K + 1)) matrix and beq as 24 X 1 matrix. For example the first row in Aeq corresponds to the 1st hour and contains the respective coefficients for the first 24 columns of the matrix followed by zeros. The second row corresponding to the 2nd hour, is basically the first row but shifted down by 24 columns. That is the first 24 columns are zeros, the next 24 columns are now the first 24 columns from the first row, followed by zeros. Also beq now contains the load demand over a 24 hour period.

Similarly, A is (num of gen) X (24 * num of gen * (K + 1)) matrix which I created using your suggestion from post on the 23 May 2018. And b is a (num of gen) X 1 matrix containing zeros.

I was able to run the simulation for this and got the optimal solutions for each hour not considering minimum up and down times (as I am yet to figure out how I should incorporate it into the code).

I have attached the code for the above if you want to have a look.

[Micah Mungal]

I was able to get the code working using the optimproblem approach. I do not think I will continue with this method. Thank you very much for your help so far.

[Alan Valadez]

Hey Micah, I would like to know how did you solve the problem, please could you show you code, i have a similar problem and I have no idea from where to start