Bi-level Optimization Problem

Question

Bashar 2023-4-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1946778-bi-level-optimization-problem

评论： Bashar 2023-4-18

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 ?

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2023-4-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1946778-bi-level-optimization-problem#answer_1215363

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.

18 个评论
显示 16更早的评论隐藏 16更早的评论

Bashar 2023-4-13

编辑：Bashar 2023-4-13

在 MATLAB Online 中打开

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 2023-4-13

编辑：Torsten 2023-4-13

在 MATLAB Online 中打开

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 

Matt J 2023-4-13

编辑：Matt J 2023-4-13

在 MATLAB Online 中打开

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

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 2023-4-15

编辑：Torsten 2023-4-15

在 MATLAB Online 中打开

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

Torsten 2023-4-17

编辑：Torsten 2023-4-17

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.

Torsten 2023-4-18

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 2023-4-18