MILP problem which is a function of time

i have solved my MILP problem using intlinprog, that part was fine. I am now extending the same problem but now it is depedent on time (i.e. I was to solve the optimization problem for each hour time step).

my question is: How do I extend my objective function and constraints to be time depedent?

thank you

[Ameer Hamza]

Which parameters of MILP change after each hour?

[sibabalo noludwwe]

equality and inequality constraints as well as the objective function. am not sure if by parameters is that what you meant

[Ameer Hamza]

If the constraints change, then you can use a for loop to solve the optimization problem after each hour.

[sibabalo noludwwe]

Wont that option have computational expense

[Ameer Hamza]

Since your constraints change, you are solving a new optimization problem. You may get some benefit by using the solution from the previous step as the initial guess for the next step.

[Ameer Hamza]

sibabalo's comment moved here to keep the discussion organized.

Maybe i did not explain it properly. The constraints equations don't change. To be exact what changes are power values from PV,wind turbines and battery storage that used in the constraints.

[Ameer Hamza]

MATLAB intlinprog function does not provide any option to deal with varying parameters. The only way I can think of is to feed the solution of the previous step as the initial guess. It will help the optimizer to converge quickly.

[sibabalo noludwwe]

Thanks ill try that. May you refer me to a matlab site of how i can achieve to use my previous results for the next if possible

[Ameer Hamza]

Please see my answer below for explanation.

Answer 1

[Ameer Hamza]

I have modified this example from the documentation: https://www.mathworks.com/help/optim/ug/intlinprog.html#mw_cd504ff9-629b-402f-acb0-7aa7c7926290 to explain the effect of using the solution from the previous time step as an initial point for the next iteration. In the following example, I used a for loop to solve an MILP problem 20 times, each time slightly varying the objective function.

First, check the code which does not use the solution from the last iteration as an initial guess of next iteration. The initial guess is always the same.

Aeq = [22    13    26    33    21     3    14    26
    39    16    22    28    26    30    23    24
    18    14    29    27    30    38    26    26
    41    26    28    36    18    38    16    26];
beq = [ 7872
       10466
       11322
       12058];
   
N = 8;
lb = zeros(N,1);
intcon = 1:N;
f = [2    10    13    17     7     5     7     3];
x0 = [8 62 23 103 53 84 46 34];
opts = optimoptions('intlinprog', 'Display', 'off');

tic
for i=1:20
    f(3:5) = randi([1 15], 1, 3); % objective function is changed in each iteration
    [x2,fval2] = intlinprog(f,intcon,[],[],Aeq,beq,lb,[],x0,opts);
end
toc

Time of execution

Elapsed time is 24.518613 seconds.

The following code update the initial guess

Aeq = [22    13    26    33    21     3    14    26
    39    16    22    28    26    30    23    24
    18    14    29    27    30    38    26    26
    41    26    28    36    18    38    16    26];
beq = [ 7872
       10466
       11322
       12058];
   
N = 8;
lb = zeros(N,1);
intcon = 1:N;
f = [2    10    13    17     7     5     7     3];
x0 = [8 62 23 103 53 84 46 34];
opts = optimoptions('intlinprog', 'Display', 'off');

tic
for i=1:20
    f(3:5) = randi([1 15], 1, 3);
    [x2,fval2] = intlinprog(f,intcon,[],[],Aeq,beq,lb,[],x0,opts);
    x0 = x2;
end
toc

Time of execution

Elapsed time is 8.489186 seconds.

This example shows a speed gain of about 3 times. The actual gain can vary from problem to problem, but still, it is better than just using a random initial guess.

[sibabalo noludwwe]

this line for example is in AMPL code: var PL {i in Lines,p in Phases, t in Time} >=-4000, <= 4000; and the constraint is time depedent and this is what I want to do in matlab for every constraint.

[Ameer Hamza]

I don't know AMPL. Can you show the mathematical form of the equation?

[sibabalo noludwwe]

please find attached

[Ameer Hamza]

Ok. So you are summing the objective function at each time step. This makes it a single optimization problem. I guess the optimization variables are , while are the weights. To simplify the notation, consider following matrices

W = [w_1 w_2 --- w_N]'; % column matrix
Gamma = [gamma_11 gamma_21 --- gamma_N1;
         gamma_12 gamma_22 --- gamma_N2;
         | | | | | | | | | | | | | | |;
         gamma_1T gamma_2T --- gamma_NT;

Using these matrices, objective function can be written as

Obj_fun = sum(Gamma*W)

After some manipulation, you can find that it is equivalent to

Obj_fun = sum(Gamma,1)*W

And this is the form of the objective function which intlinprog accepts. You can use

f = sum(Gamma,1)

as first input of intlinprog.

Similarly, you can write down all the constraints and create the matrix Aeq and Beq using the cofficients of .

[sibabalo noludwwe]

I have done this for the part where the optimization problem is not time dependent. what I need now is icoporating a time index on my variables and that is what am struggling with and have been trying to exlpain.

[sibabalo noludwwe]

can I still use the same for loop if i am using CPLEX to solve the MILP problem

[Ameer Hamza]

I haven't use CPLEX, but I just had a brief look. If it allows a similar interface, then you should be able to use it in place of intlinprog.

[sibabalo noludwwe]

I think It does.

With the initial solution. Am I right if I run the simulation outside the for loop and then use that as my x0.?

[Ameer Hamza]

Yes, that should also work. In fact it can be a good thing since it can give a good starting point to the optimization algorithm.

[sibabalo noludwwe]

is it posible that i am doing something wrong, my objective function is same for all the the time steps but when am running the simulations seperately they are not the same.

[Ameer Hamza]

The output of an optimizer can depend on its initial guess. Maybe that is why you are getting different solutions.

[sibabalo noludwwe]

thanks

[sibabalo noludwwe]

am now experiencing an odd problem when i am running the optimisation problem for every time step seperately I get an optimal solution but when I use a for loop for some reason there is no optimal solution even when I dont add the initial solution.

[Ameer Hamza]

Can you add a breakpoint in your code with for-loop and run it step by step using the debugging options available in the toolbar. It will help you to find out which line is not working as expected.

[Ameer Hamza]

sibabalo's comment posted as answer moved here:

should I add the breakpoint inside the for loop. coz whats happening is that it runs for the first round and gives an optimal solution and stills runs all the times steps but no with no optimal solution. so am not sure if that explaination is for this. because the for loops runs but the intlinprog returns no optimal solution when I run the for loop.

[Ameer Hamza]

Are you changing the objective function of MILP in each iteration or solving the same problem in each iteration.

[sibabalo noludwwe]

Am solving the same objective function the only thing that changes is the ub.

for i=1:10

ubpv1=[zeros(33-29,1);pvpower(S1,2);zeros(33-5,1)];

ubpv2=[zeros(33-2,1);pvpower(S1,2)+windpower(S1,2);zeros(33-32,1)]; %%

ub=[ub1;ubp1;ubp2;ubq1;ubq2;ubdgp1;ubdgp2;ubdgq1;ubdgq2;ubpv1;ubpv2;ubv12;ubpm];

this is what is changing only for every S1 value.

[x,fval,EXITFLAG,OUTPUT] = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub,options);

fprintf ('Number of critical loads restored = %f\n',-1*fval);

disp ('Control variables =');

disp (x');

S1=S1+1;

end

my confusion is that when i change the value of S1 manually and remove the for loop I get optimal solution...

[Ameer Hamza]

This is quite strange behavior. I couldn't see any reason why this should be happening. I recommend starting two instances of MATLAB and run both code line by line. You will definitely see the difference at some point.

[sibabalo noludwwe]

I am running the same code i just comment the for loop and bring it back in. So its the code

Answer 2

[sibabalo noludwwe]

LP: Optimal objective value is -4.403333.

Branch and Bound:

nodes total num int integer relative

explored time (s) solution fval gap (%)

3 1.95 1 -3.750000e+00 9.403509e+00

4 2.67 2 -3.930000e+00 5.409060e+00

5 3.13 2 -3.930000e+00 0.000000e+00

Optimal solution found.

Intlinprog stopped because the objective value is within a gap tolerance of the optimal value,

options.AbsoluteGapTolerance = 0 (the default value). The intcon variables are integer within tolerance,

options.IntegerTolerance = 1e-05 (the default value).

Number of critical loads restored = 3.930000

Control variables =

Columns 1 through 11...... and so on

for the other time ste this is the solution

No feasible solution found.

Intlinprog stopped because no point satisfies the constraints.

Number of critical loads restored =

Control variables =