Bi-level Optimization Problem

Hello all

I have to solve a bi-level optimization problem that is described as follows:

The upper layer problem is needed to be solved by using 'ga' solver.

and the lower problem is a MILP that should be solved by 'intlinprog'

The procedure of the solving the problem is as follows (in general):

1st ) the upper decision variables should be initialized and passed to lower problem as fixed values to start solving it .

2nd ) the lower problem is run and solved and produce the solutions ( the lower decision variables)

3rd ) the lower problem solutions should be returned to upper problem ( as values ) to start solving the upper problem and producing the solutions.

4th) the resultant upper solutions should be paassed to lower problem again until the maximum number of ga's generation is reached

The figures attached below clarify the structure of the mentioned problem

These figures are taken from this paper (https://www.sciencedirect.com/science/article/abs/pii/S0960148119303799)

My questions are 1) how can I pass the intial populations (initial values of upper decision variables) to start the lower porblem ?

2)Also, how can I pass the solution of the upper problem to lower problem after each generation of ga ?

Answer 1

[Matt J]

Sine intlinprog is called by by ga's fitness function, you will have the variables in the upper problem as input to the lower problem there.

[Bashar]

Thank you for your fast response @Matt J

Let me clarify more, the variables in upper problem should be initialized ( the initial populations) then passed to lower problem as fixed values. The lower (inner loop) problem has own variables and totally separated from those of upper problem. Afterthat, the lower problem is run by intlnprog and produces the solutions ( values of lower variables ), Here, I know how to pass them to upper problem. but, my question how can I pass the upper variables as fixed values in these cases

1) before starting the lower problem (passing the initial populations)

2) and after each generation of ga ( passing the solutions to the lower problem again)

in other words, which option related to 'Ga' toolbox can achieve this task ?

I hope I could clarify my point and didn't confuse your understanding the question

Thanks in advance

[Torsten]

You don't need to care about passing variables or something like this.

The interface where upper and lower problem meet is the function "nonlcon" called by the upper problem where you define its constraints. Here, the solution variables of the upper problem so far are available, and you have to use them to solve the lower problem. With this solution, you can define the constraints of the upper problem.

If it's additionally necessary to get the solution of the lower problem to define the objective of the upper problem, you have to do the same in the function where you define this objective for ga. If this is the case, you might be able to save optimizations with intlinprog following the procedure described here:

https://de.mathworks.com/help/optim/ug/objective-and-constraints-using-common-function.html

Maybe it's possible also to use the problem-based approach to solve this problem, but I'd prefer the solver-based approach in this complicated case.

[Bashar]

Thank you @Torsten for your valued response

Carefully, I read your comment and procedure in the link and I tried to understand what you meant exactly, however, I got confused, I'm so sorry .

Also, I'm familiar with problem-based approach, I wrote the script of the lower problem alone. I have never tried the solver-based approach

but below I povide an example (script) I hope that clarifies my problem

the main script

clc 
clear 
close 
t = 1; % represent a fixed parameter. 
[UpperSol,Upperfval,exitflag_Upper,LowerSol,Lowerfval,exitflag_Lower] = Test2(t)

and main code ( lower and upper problem ) is presented here

function [UpperSol,Upperfval,exitflag_Upper,LowerSol,Lowerfval,exitflag_Lower] = Test2(t) 
x = optimvar("x","LowerBound",0,"UpperBound",1,"Type","integer");
y = optimvar("y","LowerBound",0,"UpperBound",20,"Type","integer");
f = optimvar("f");
camxy = @(x,y,f)(4 - 2.1.*x.^2 + x.^4./3).*x.^2 + x.*y + (-4 + 4.*y.^2).*y.^2 + f;
cam = camxy(x,y,f);
[LowerSol,Lowerfval,exitflag_Lower] = nestfun1(x,y);

 function [LowerSol,Lowerfval,exitflag_Lower] = nestfun1(x,y) 
 a = optimvar('a','LowerBound',0);
 b = optimvar('b','LowerBound',0);

c = optimvar('c','LowerBound',0,'UpperBound',1,'Type','integer');
lowerProb = optimproblem('Objective',-3*a*x-2*b-c);

lowerProb.Constraints.cons1 = a*y + b + c <= 7;
lowerProb.Constraints.cons2 = 4*a + 2*b + c == 12;

    [lowerSol,Lowerfval,exitflag_Lower] = solve(lowerProb)
    LowerSol = struct2cell(lowerSol);
   end 
upperProb = optimproblem("Objective",cam);
upperProb.Constraints.cons1 = t*x*y + x - y + 1.5 <= 0;
upperProb.Constraints.cons2 = 10 - x*y <= 0;
upperProb.Constraints.cons3 = f == Lowerfval;
options = optimoptions(@ga,...
    'PlotFcn',{@gaplotbestf,@gaplotmaxconstr},...
    'Display','iter');
[upperSol,Upperfval,exitflag_Upper] = solve(upperProb,"Solver","ga","Options",options);
UpperSol = struct2cell(upperSol);
end

in script above, the (x and y) are upper problem's variables

and the (a, b and c) represent the lower problem's variables .

when the x and y paased to lower problem as fixed values ( x and y don't change their values during the lower problem ), the lower problem becomes as MILP problem ( linear objective function and constraints included integer variable (c) )

kindly, explain your answer by this example, please.

[Matt J]

Test2(1)

Solving problem using intlinprog.
LP:                Optimal objective value is -12.000000.                                           


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05.

Error using constrValidate
GA does not solve problems with integer and equality constraints.
For more help see No Equality Constraints in the documentation.

Error in gacommon (line 185)
[ConstrInfo,Iterate,output] = constrValidate(NonconFcn, ...

Error in ga (line 377)
NonconFcn,options,Iterate,type] = gacommon(nvars,fun,Aineq,bineq,Aeq,beq,lb,ub, ...

Error in solution>Test2 (line 10)
ga(@cam, 3, [],[],[],[],[0 0 -inf], [1,20,+inf],@(p)nonlcon(p,t),1:2,options);

function [UpperSol,Upperfval,exitflag_Upper,LowerSol] = Test2(t) 

  options = optimoptions(@ga,...
    'PlotFcn',{@gaplotbestf,@gaplotmaxconstr},...
    'Display','iter');

     [p, Upperfval,exitflag_Upper] = ...
         ga(@cam, 3, [],[],[],[],[0 0 -inf], [1,20,+inf],@(p)nonlcon(p,t),1:2,options);

      UpperSol.x=p(1);
      UpperSol.y=p(2);
      UpperSol.f=p(3);
     [~,~,LowerSol,Lowerfval,exitflag_Lower] = nonlcon(UpperSol,t);

end

function out=cam(p)

    [x,y,f]=deal(p(1),p(2),p(3));
    
    out=(4 - 2.1.*x.^2 + x.^4./3).*x.^2 + x.*y + (-4 + 4.*y.^2).*y.^2 + f


end

function [cineq,ceq,LowerSol,Lowerfval,exitflag_Lower] = nonlcon(p,t)

    [x,y,f]=deal(p(1),p(2),p(3));

           %%%Sub-problem (START)
             a = optimvar('a','LowerBound',0);
             b = optimvar('b','LowerBound',0);  
             c = optimvar('c','LowerBound',0,'UpperBound',1,'Type','integer');
     
            lowerProb = optimproblem('Objective',-3*a*x-2*b-c);    
            lowerProb.Constraints.cons1 = a*y + b + c <= 7;
            lowerProb.Constraints.cons2 = 4*a + 2*b + c == 12;
            
             [LowerSol,Lowerfval,exitflag_Lower] = solve(lowerProb);
             
           %%%Sub-problem (END)
           
    cineq(1) = t*x*y + x - y + 1.5 ;
    cineq(2) = 10 - x*y ;  
    ceq= f - Lowerfval ;
   
     if exitflag_Lower<0
     
       ceq=inf;
       
     end
 
end 

[Torsten]

ga restriction makes it impossible to solve the test problem.

[Bashar]

Thanks for your answers @Matt J and @Torsten

Actually, I combined lower and upper problems in this example from different sources form mathworks' decumentations, in order to formulate an example clairfied my question; so, you can find both examples from these sources MILP Example and Ga (upper problem) example

So, I think the problem is the combination of these two examples is inappropriate (in my opinion )

Anyway, I will try your method in my real problem and I will give you the feedback here

Thank you so much for your efforts

Kind Regards

[Matt J]

Here's an alternative formulation of the example problem that avoids the previously encountered errors, though for t=1 a solution may not exist.

Test2(1)

Single objective optimization:
3 Variable(s)
2 Integer variable(s)
2 Nonlinear inequality constraint(s)

Options:
CreationFcn:       @gacreationuniformint
CrossoverFcn:      @crossoverlaplace
SelectionFcn:      @selectiontournament
MutationFcn:       @mutationpower

                                  Best          Mean         Stall
Generation      Func-count     Penalty         Penalty    Generations
    1               80          0.1429          0.5096        0
    2              118          0.1429          0.2911        1
    3              156          0.1429          0.2818        2
    4              194          0.1429          0.1979        3
    5              232          0.1429          0.1643        4
    6              270          0.1429          0.2004        5
    7              308          0.1429          0.1754        6
    8              346          0.1429          0.1718        7
    9              384          0.1429          0.2218        8
   10              422          0.1429          0.1968        9
   11              460          0.1429          0.1693       10
   12              498          0.1429          0.1529       11
   13              536          0.1429           0.185       12
   14              574          0.1429          0.1893       13
   15              612          0.1429          0.2054       14
   16              650          0.1429          0.1868       15
   17              688          0.1429          0.2179       16
   18              726          0.1429          0.2118       17
   19              764          0.1429          0.2207       18
   20              802          0.1429          0.2243       19
   21              840          0.1429          0.2296       20
   22              878          0.1429          0.1943       21
   23              916          0.1429          0.1729       22
   24              954          0.1429          0.1729       23
   25              992          0.1429          0.1929       24
   26             1030          0.1429          0.2107       25
   27             1068          0.1429          0.1654       26
   28             1106          0.1429          0.1818       27
   29             1144          0.1429          0.2218       28

                                  Best          Mean         Stall
Generation      Func-count     Penalty         Penalty    Generations
   30             1182          0.1429          0.1654       29
   31             1220          0.1429          0.1454       30
   32             1258          0.1429          0.1704       31
   33             1296          0.1429          0.1889       32
   34             1334          0.1429          0.1479       33
   35             1372          0.1429          0.1604       34
   36             1410          0.1429          0.1604       35
   37             1448          0.1429          0.1918       36
   38             1486          0.1429          0.1554       37
   39             1524          0.1429          0.1454       38
   40             1562          0.1429          0.2218       39
   41             1600          0.1429          0.1954       40
   42             1638          0.1429          0.2193       41
   43             1676          0.1429          0.1629       42
   44             1714          0.1429          0.2257       43
   45             1752          0.1429          0.1793       44
   46             1790          0.1429          0.1629       45
   47             1828          0.1429          0.1893       46
   48             1866          0.1429          0.1904       47
   49             1904          0.1429          0.1957       48
   50             1942          0.1429          0.2207       49
   51             1980          0.1429          0.1843       50

Solving problem using intlinprog.
LP:                Optimal objective value is -12.000000.                                           


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05.

Optimization terminated: average change in the penalty fitness value less than options.FunctionTolerance
but constraints are not satisfied.

Solving problem using intlinprog.
LP:                Optimal objective value is -12.000000.                                           


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05.

ans = struct with fields:
    x: 1
    y: 16
    f: -9.3256e+03

function [UpperSol,Upperfval,exitflag_Upper,LowerSol] = Test2(t) 

  options = optimoptions(@ga,...
    'PlotFcn',{@gaplotbestf,@gaplotmaxconstr},...
    'Display','iter');

     [p, Upperfval,exitflag_Upper] = ...
         ga(@cam, 3, [],[],[],[],[0 0 -inf], [1,20,+inf],@(p)nonlcon(p,t),1:2,options);

      UpperSol.x=p(1);
      UpperSol.y=p(2);
      UpperSol.f=p(3);
     [~,LowerSol,Lowerfval,exitflag_Lower] = cam(p);

end

function [out,LowerSol,Lowerfval,exitflag_Lower]=cam(p)

    [x,y,f]=deal(p(1),p(2),p(3));
    
             %%%Sub-problem (START)
             a = optimvar('a','LowerBound',0);
             b = optimvar('b','LowerBound',0);  
             c = optimvar('c','LowerBound',0,'UpperBound',1,'Type','integer');
     
            lowerProb = optimproblem('Objective',-3*a*x-2*b-c);    
            lowerProb.Constraints.cons1 = a*y + b + c <= 7;
            lowerProb.Constraints.cons2 = 4*a + 2*b + c == 12;
            
             [LowerSol,Lowerfval,exitflag_Lower] = solve(lowerProb);
             
           %%%Sub-problem (END)
    
           
    if  exitflag_Lower>=0
     out=(4 - 2.1.*x.^2 + x.^4./3).*x.^2 + x.*y + (-4 + 4.*y.^2).*y.^2 + Lowerfval;
    else
     out=inf;   
    end

end

function [cineq,ceq] = nonlcon(p,t)

    [x,y,f]=deal(p(1),p(2),p(3));

           
    cineq(1) = t*x*y + x - y + 1.5 ;
    cineq(2) = 10 - x*y ;  
    ceq=[];
 
end 

[Torsten]

@Bashar

So, I think the problem is the combination of these two examples is inappropriate (in my opinion )

No. The problem was that "ga" does not accept nonlinear equality constraints if the problem also has integer constraints:

https://de.mathworks.com/help/gads/ga.html

[Bashar]

@Matt J and @Torsten

Excuse me, I have question

if I want to use the lower solutions (a,b and c as values) in adding third ineqaulity constraint for the upper problem as follow

cineq(3) = x*a + b*y;

How can I modify the script above ?

[Torsten]

There might be potential to save time in exchanging solutions of your lower problem between "cam" and "nonlcon" with the help of the link I already gave you:

https://de.mathworks.com/help/optim/ug/objective-and-constraints-using-common-function.html

Adding the constraint to the lower problem needs less changes to the code from above than adding it to the upper problem (as you want).

I hope you are aware that with this constraint, your problem becomes infeasible.

Test2(1)

function [UpperSol,Upperfval,exitflag_Upper,LowerSol] = Test2(t) 

  options = optimoptions(@ga,...
    'PlotFcn',{@gaplotbestf,@gaplotmaxconstr},...
    'Display','iter');

     [p, Upperfval,exitflag_Upper] = ...
         ga(@cam, 3, [],[],[],[],[0 0 -inf], [1,20,+inf],@(p)nonlcon(p,t),1:2,options);

      UpperSol.x=p(1);
      UpperSol.y=p(2);
      UpperSol.f=p(3);
     [~,LowerSol,Lowerfval,exitflag_Lower] = cam(p);

end

function [out,LowerSol,Lowerfval,exitflag_Lower]=cam(p)

    [x,y,f]=deal(p(1),p(2),p(3));
    
             %%%Sub-problem (START)
             a = optimvar('a','LowerBound',0);
             b = optimvar('b','LowerBound',0);  
             c = optimvar('c','LowerBound',0,'UpperBound',1,'Type','integer');
     
            lowerProb = optimproblem('Objective',-3*a*x-2*b-c);    
            lowerProb.Constraints.cons1 = a*y + b + c <= 7;
            lowerProb.Constraints.cons2 = 4*a + 2*b + c == 12;
            
             [LowerSol,Lowerfval,exitflag_Lower] = solve(lowerProb);
             
           %%%Sub-problem (END)
    
           
    if  exitflag_Lower>=0
     out=(4 - 2.1.*x.^2 + x.^4./3).*x.^2 + x.*y + (-4 + 4.*y.^2).*y.^2 + Lowerfval;
    else
     out=inf;   
    end

end

function [cineq,ceq] = nonlcon(p,t)

    [x,y,f]=deal(p(1),p(2),p(3));

        [x,y,f]=deal(p(1),p(2),p(3));
    
             %%%Sub-problem (START)
             a = optimvar('a','LowerBound',0);
             b = optimvar('b','LowerBound',0);  
             c = optimvar('c','LowerBound',0,'UpperBound',1,'Type','integer');
     
            lowerProb = optimproblem('Objective',-3*a*x-2*b-c);    
            lowerProb.Constraints.cons1 = a*y + b + c <= 7;
            lowerProb.Constraints.cons2 = 4*a + 2*b + c == 12;
            
             [LowerSol,Lowerfval,exitflag_Lower] = solve(lowerProb);
             
           %%%Sub-problem (END)
           
    cineq(1) = t*x*y + x - y + 1.5 ;
    cineq(2) = 10 - x*y ;  
    cineq(3) = LowerSol.a*x + LowerSol.b*y;
    ceq=[];
 
end

[Bashar]

Thanks a lot @Torsten

Your response is always valued .

I thought in this way of modification, solving the lower problem twice, but i didn't know how

I have read the documentation inside the link, I think that using the Parallel Computing Toolbox will reduce the execution time. I have run the code with and without this option, the result was great, the run time were declined by around 40 sec .

With respect to this constraint, I added it arbitrary in order to simulate my actual problem (you can say as a test ) because my problem is larger w.r.t the number of constraints ,varaibles and input data .

I appreciate your help and Mr @Matt J

Kind regards

[Matt J]

@Bashar is your problem solved, then?

[Bashar]

@Matt J

If you mean My actual problem, currently I'm working on it as your and @Torsten's instructions in this test problem.

Thanks a lot

[Torsten]

Hello @Matt J and @Torsten

Kindly, Could anyone help me with this hypothetical case in the above test code ?

the case is the following: If I want to add a forth variable (z) (vector (state) variable with lenght of N) to the upper problem, and this variable is only incorporated in the upper problem as (state) updating constraint as follows

z(t+1) = z(t) + a(t)/b; for t=1:N-1 and a and b are the lower problem variables ( as you know)

this is required that the lower variable (a) should be redfined in lower problem as a vector with length of N

With regardless of the modification in the lower problem, I need just how can I define the variable (z) in the problem, especially in nonlcon function where the constraints are defined.

I searched in Mathworks documentation for longtime and I didn't find anything may help me; in particular, how to define vector variable in the solver based appraoch.

Actually, I introduce this case to simulate challenges will appear in my actual problem ( which I can't share it here).

I appreciate your valued help

Kind regards

[Torsten]

Everything about the direct call of ga with many examples is explained here:

https://de.mathworks.com/help/gads/ga.html

At the moment, you call ga as

[p, Upperfval,exitflag_Upper] = ...

ga(@cam, 3, [],[],[],[],[0 0 -inf], [1,20,+inf],@(p)nonlcon(p,t),1:2,options);

So you could change the 3 to N+3,

the upper and lower limits to [0 0 -inf zlb] resp. [1,20,+inf,zub] where zlb and zub are row vectors of length N which contain the lower and upper bounds for the z components

and maybe

1:2 to a modified index vector with the possible z components that are integer variables.

But note that this way will not work because you cannot define z as optimization variables and the relation for z in the upper problem because you had to define the relations as nonlinear equality constraints. But as you might remember from your first question here, ga does not accept nonlinear equality constraints together with integer constraints.

I suggest you don't stick to your test problem too strong because it only shows you how in principle a coupling of two optimizations works. You should now concentrate on your real problem and plan which variables to solve in the upper, which in the lower problem, which constraints are to be set in the upper and which in the lower problem and if this can be done with the MATLAB solvers you plan to use.

[Bashar]

According to your last paragraph, I totally agree with you. Currently, I'm working on my real problem, However, I always would like to predict the challanges and difficulties before starting any important work. This is my way in dealing with such tasks.

I’m sorry for asking so many questions. I appreciate your time in the responding . There are no words to express my thanks.

Kind Regards

[Torsten]

I'm also a fan of "work from easy to difficult". But in your case, all the preparatory work has been done, you know how to couple two optimizations and it doesn't make sense to extend the test problem. You will have to start anew with your real problem making use of the method, but not the equations of the test problem.

[Bashar]

You're completely right. Absolutly, I shall utilze the method not the same equations of the test problem.

Mnay Thanks