Solving Multi integer problem with non linear constraints

Hello, I am trying to solve a multi integer programming problem that has some non linear constraints in it. The problem basically comprises of 7 continuous variables that are evaluated separately on a 24 hour interval, hence an altogether of 168 variables (7x24=168), and 4 binary variables which also are to be evaluated on the same interval (24x7= 96 integer variables). 1) First I tried solving it with genetic algorithm (ga) solver, but it was giving error for the non linear equality constraints. So I transformed each non linear equality constraint to 2 inequality constraints that bounds the equality constraint in a certain tolerance. The problem was that the final solution was not respecting the equality constraints. 2) Then I tried using surrogateopt solver, which apparently is designed to handle the non linear equality constraints, but even the solution obtained from it did not respecte the equality constraints. 3) Finally I thought of using some pen and paper, and noticed that there were 5 equality constraints with 7 continuous variables in it. I transformered 5 of the continuous variables as a function of 2 continous variables, and placed the corresponding expression for each unknown variable using optimexpr(). This procedure allowed me to decrease the number of variables from a total of 264 (168+96) variables to 144 variables .Naturally to place bounds for those 5 continuous variables I needed to use inequality constraints. Using this approach the equality constraints seemed to work, but now the problem was that none of the inequality constraint was being respected by the either of the aforementioned solvers. For genetic algorithm I increased the population size upto 7000 along with max generations of 1000 ,which took over 7 hours to solve, and for surrogate I used MaxFunctionEvaluations of over 1000, which eventually stopped working at about 700 evaluations. Now my question is that is MATLAB’s Optimization toolbox capable to solve multi integer non linear problem with about 150 variables? If yes, what sort of modifications can I bring in my code, if there are not any, is there another solver other than the ones I used that can help me? Because in one of the post on this forum I read that both the solvers works best with 100 variables, though it also mentioned there is no defined limit. P.S: I used problem based approach for each of the solver. Thanks in advance! Basil

[Matt J]

Kanwar Basil Babar's comment moved here:

clc;clearvars

ga_trial=optimproblem('ObjectiveSense','minimize');
%DATA
Ql=[0,0,0,0,0,10,10,8,6,4,4,4,3,3,4,6,8,10,10,10,8,4,0,0]';
Tgb_s=100;
m_L_max=0.2;
m_gb_min=0.05*1;
m_gb_max=0.3*1;
Tr_min=40;
dT=1;
Ms=m_gb_max*dT; 
Mr=Ms;
As= 5*(Ms/1000)^(2/3); 
Qgb_max=24;
Qgb_min=6;
c=4.178;
T_amb=25;
h=0.040/0.050/1000;
eta_gb=0.9;
f_ts=2/60;
tol=1e-3;

m_L=zeros(24,1); %Vector for m_L in operational hours

for i=1:1:24
    if Ql(i)>0
        m_L(i)=m_L_max;
    else
        m_L(i)=0;
    end
end

%Variables

%Continuous Variables
Qgb= optimvar('Qgb',24,'Type','continuous','LowerBound',0,'UpperBound',Qgb_max); 
m_gb=optimvar('m_gb',24,'Type','continuous','LowerBound',0,'UpperBound',m_gb_max);
m_sr=optimvar('m_sr',24,'Type','continuous','LowerBound',0,'UpperBound',m_gb_max);
m_rs=optimvar('m_rs',24,'Type','continuous','LowerBound',0,'UpperBound',m_L_max);
Ts=optimvar('Ts',24,'Type','continuous','LowerBound',0,'UpperBound',Tgb_s);
Tr=optimvar('Tr',24,'Type','continuous','LowerBound',Tr_min,'UpperBound',Tgb_s); 
TL_r=optimvar('TL_r',24,'Type','continuous','LowerBound',0,'UpperBound',Tgb_s);

%Integer Variables
delta_rs=optimvar('delta_rs',24,'Type','integer','LowerBound',0,'UpperBound',1);
delta_gb=optimvar('delta_gb',24,'Type','integer','LowerBound',0,'UpperBound',1); 
z=optimvar('z',24,'Type','integer','LowerBound',0,'UpperBound',1);
delta_sr=optimvar('delta_sr',24,'Type','integer','LowerBound',0,'UpperBound',1);

%CONSTRAINTS

%Energy Balance for Boiler 
con1a=optimconstr(24,1); %Non Linear Constraint
con1b=optimconstr(24,1); %Non Linear Constraint
for i=1:1:24
    con1a(i)=Qgb(i)<=m_gb(i)*c*(Tgb_s-Tr(i))+tol;
    con1b(i)=Qgb(i)>=m_gb(i)*c*(Tgb_s-Tr(i))-tol;
end

ga_trial.Constraints.con1a = con1a;
ga_trial.Constraints.con1b = con1b;


%Operational Limit for the boiler
con2a= delta_gb*Qgb_min<=Qgb;
con2b= delta_gb*Qgb_max>=Qgb;

%Penalty Function for the Boiler
con2c= delta_gb(2:24)-delta_gb(1:23)<=z(2:24);
con2d= delta_gb(1)-delta_gb(24)<=z(1);


ga_trial.Constraints.con2a = con2a;
ga_trial.Constraints.con2b = con2b;
ga_trial.Constraints.con2d = con2c;
ga_trial.Constraints.con2d = con2d;

%Operational Limit for the mass flow through the boiler
con3a= delta_gb*m_gb_min<=m_gb;
con3b= delta_gb*m_gb_max>=m_gb;

ga_trial.Constraints.con3a = con3a;
ga_trial.Constraints.con3b = con3b;

%Energy Balance for the Load
con4=optimconstr(24,1);
for i=1:1:24
    con4(i)= Ql(i)-m_L(i)*c*(Ts(i)-TL_r(i))==0;
end
ga_trial.Constraints.con4 = con4;

%Operational Limit for mass flow from Tank-S to Tank-R && Tank-R to Tank-S
con5a= m_sr<=delta_sr*m_gb_max;
con5b= m_rs<=delta_rs*m_L_max;
con5c= delta_rs+delta_sr<=1;

ga_trial.Constraints.con5a = con5a;
ga_trial.Constraints.con5b = con5b;
ga_trial.Constraints.con5c = con5c;

%Mass Balance on either of the Tank
con6= m_sr+m_L-m_rs-m_gb==0;

ga_trial.Constraints.con6 = con6;

%Energy Balance on Tank-S
con7a_empty=optimexpr(23,1); %Non Linear Constraint
for i=1:1:23
    con7a_empty(i)= Ts(i+1)-(1/(Ms*c)*(Ms*c*Ts(i)+m_gb(i)*c*Tgb_s*dT+m_rs(i)*c*Tr(i)*dT-m_L(i)*c*Ts(i)*dT-m_sr(i)*c*Ts(i)*dT-h*As*(Ts(i)-T_amb)*dT));
end
con7a1= con7a_empty>=tol;
con7a2= con7a_empty<=-tol;
con7b1= Ts(1)>=1/(Ms*c)*(Ms*c*Ts(24)+m_gb(24)*c*Tgb_s*dT+m_rs(24)*c*Tr(24)*dT-m_L(24)*c*Ts(1)*dT-m_sr(24)*c*Ts(24)*dT-h*As*(Ts(24)-T_amb)*dT);
con7b2= Ts(1)<=1/(Ms*c)*(Ms*c*Ts(24)+m_gb(24)*c*Tgb_s*dT+m_rs(24)*c*Tr(24)*dT-m_L(24)*c*Ts(1)*dT-m_sr(24)*c*Ts(24)*dT-h*As*(Ts(24)-T_amb)*dT);

ga_trial.Constraints.con7a1 = con7a1;
ga_trial.Constraints.con7a2 = con7a2;
ga_trial.Constraints.con7b1 = con7b1;
ga_trial.Constraints.con7b2 = con7b2;

%Energy Balance on Tank-R
con8a_empty=optimexpr(23,1); %Non Linear Constraint
for i=1:1:23
    con8a_empty(i)= Tr(i+1)-(1/(Mr*c)*(Mr*c*Tr(i)+m_sr(i)*c*Ts(i)*dT+m_L(i)*c*TL_r(i)*dT-m_gb(i)*c*Tr(i)*dT-m_rs(i)*c*Tr(i)*dT-h*As*(Tr(i)-T_amb)*dT));
end
con8a1= con8a_empty>=0;
con8a2= con8a_empty<=0;
con8b1= Tr(1)>=1/(Mr*c)*(Mr*c*Tr(24)+m_sr(24)*c*Ts(24)*dT+m_L(24)*c*TL_r(24)*dT-m_gb(24)*c*Tr(24)*dT-m_rs(24)*c*Tr(24)*dT-h*As*(Tr(24)-T_amb)*dT);
con8b2= Tr(1)<=1/(Mr*c)*(Mr*c*Tr(24)+m_sr(24)*c*Ts(24)*dT+m_L(24)*c*TL_r(24)*dT-m_gb(24)*c*Tr(24)*dT-m_rs(24)*c*Tr(24)*dT-h*As*(Tr(24)-T_amb)*dT);

ga_trial.Constraints.con8a1 = con8a1;
ga_trial.Constraints.con8a2 = con8a2;
ga_trial.Constraints.con8b1 = con8b1;
ga_trial.Constraints.con8b2 = con8b2;


%Objective Function

Fgb=Qgb/eta_gb + z*Qgb_max*f_ts/eta_gb;
obj_func=sum(Fgb);
ga_trial.Objective = obj_func;

show(ga_trial)

%Solver

opts = optimoptions(@ga, ...
                    'PopulationSize', 500, ...
                    'MaxGenerations', 5, ...
                    'EliteCount', 10, ...
                    'FunctionTolerance', 1e-8, ...
                    'PlotFcn', @gaplotbestf);


rng default % For reproducibility
[sol,fval,exitflag] = solve(ga_trial,"Solver","ga","Options",opts)

Answer 1

[Matt J]

Nothing can really be said without seeing a mathematical description of the problem. However, it doesn't seem theoretically likely that you're going to find a solution if the nonlinear equality constraints involve the binary variables. Pretty much the only time you can hope that equality-constrained integer programs will even be feasible is if they involve polynomial functions with integer coefficients, e.g.,

3*x^2+ 4*y^3=7

Also, even when this is the case, if x and y are binary variables, there's no need to be writing the constraint non-linearly. It can be replaced with a linear constraint.

3*x+4*y=7

[Kanwar Basil Babar]

I have attatched the code in my previous comment so you can get an idea about the problem, and have also identified the constraints that are non linear. I have just shown the particular code I used, without removing the non linear equality constraints.

In my case the continous variables are involved with the non linear equality constraints, so do you think there is workaround for it?

[Matt J]

Yes. For each fixed set of values of your continuous variables, you obtain an integer linear programming problem in the binary variables. In fact, it appears that the linear constraints involving the integer variables are linear in all the variables, even the continuous ones. So, you can eliminate the integer variables from the problem by solving an integer programming sub-problem inside your fitness function,

function [fval,xbinary]=fitnessFcn(xcontinuous)

  f=___ %depends on xcontinuous
  A=___    
  b=___  %depends on xcontinuous
  Aeq=___
  beq=___ %depends on xcontinuous
  lb=___ %depends on xcontinuous
  ub=___ %depends on xcontinuous
  
  [xbinary,fval,exitflag]=intlinprog(f,intcon,A,b,Aeq,beq,lb,ub);
  
  if exitflag<=0, fval=inf; end
  
end

If you are to continue to use the problem-based framwork to solve the outer optimization problem, you will need to implement the fitnessFcn() in the above functional form, and so you will need to use fcn2optimexpr(). There should be no integer variables in the outer problem, just xcontinuous.

To set-up the inner integer linear program, you might be able to use the problem-based framework as well, though you would have to generate a new optimproblem in every call to fitnessFcn(), which will be slower than just building the matrices directly.