Linear Program through (FminCon,Ga,Linprog)

I have a simple linear problem including non-linear constraints. I have tried a lot of solvers including (fmincon,ga) algorithms. But I can not achieve desired results. If I exclude nonlinear constraints only linprog can solve it to optimality. But as far as I understand fmincon and ga also can handle simple linear problems. And my question is: How to correctly configure these solvers and algorithms to solve simple linear problem and after that include nonlinear constraints to achieve feasible result. I have read User Guides for Toolboxes but I still do not understand how these solvers works. Thank you

Answer 1

[Alan Weiss]

I have no idea what your problem is. There are a lot of examples that use nonlinear constraints, including nonlinear inequality constraints, nonlinear equality and inequality constraints, and a host of other examples.

If you have a linear objective function then your objective function definition is

 fun = @(x)f'*x;

where f is a column vector of coefficients and x0 is a column vector of initial guesses. (For ga you need x0 to be a row vector and your objective can be x*f.)

Without more detail from you I can provide only generic help. For more detailed help, provide a more detailed question.

Alan Weiss

MATLAB mathematical toolbox documentation

[Anton Nester]

This is my model and output of the solvers! This is simplified model. Model:

List of variables.

        variables = {'f16','f25','f35','f47','f56','f67','f78',...
          'f810','f1012','f1213','f1216',...
          'f1617','f1718','f1819','f1921','f1923',...
          'f2324','f2326','f2628','f2830','f2832',...
          'f89','f1011','f12fd','f13fd','f1314','f1215','f16fd','f17fd',...
          'f18fd','f19fd','f21fd','f2122','f1920','f24fd','f2425','f2627','f28fd','f2829','f30fd','f3031','f32fd','f3233',...
          'delta1minus','delta2minus','delta3minus','delta4minus','delta1plus','delta2plus','delta3plus','delta4plus'};
      N = length(variables);
      % create variables for indexing
      for v = 1:N
         eval([variables{v},' = ', num2str(v),';']);
      end
Where:

Sets:

1. G-set of gas fields(In my case 1,2,3,4-four production fields) ;
2. M-set of customer nodes "In the list of var for x(f89) to x(f3233)";
3. N-set of all nodes;
4. K-set of nodes with flow and demand;
5. Fij - flow from node i to node j;
6. Fjfd - nodes with flow and demand at the same time;
7. delta"i"minus - Negative deviation from the planned production in field g;
8. delta"i"plus - Positive deviation from the planned production in field g.

Objective for "ga" solver or "fmincon";

Z=-x(f89)*40-x(f1011)*40-x(f12fd)*40-x(f13fd)*40-x(f1314)*40-x(f1215)*40-x(f16fd)*40-x(f17fd)*40-x(f18fd)*40-x(f19fd)*40-x(f1920)*40-...
   x(f21fd)*40-x(f2122)*40-x(f24fd)*40-x(f2425)*40-x(f2627)*40-x(f28fd)*40-x(f2829)*40-x(f30fd)*40-x(f3031)*40-...
    x(f32fd)*40-x(f3233)*40+x(delta1plus)*4+x(delta2plus)*4+x(delta3plus)*4+x(delta4plus)*4+x(delta1minus)*4+x(delta2minus)*4+x(delta3minus)*4+x(delta4minus)*4;

Objective for "linprog";

w([f89,f1011,f12fd,f13fd,f1314,f1215,f16fd,f17fd,f18fd,f19fd,f1920,f21fd,f2122,f24fd,f2425,f2627,f28fd,f2829,f30fd,f3031,f32fd,f3233,...
    delta1plus,delta2plus,delta3plus,delta4plus,delta1minus,delta2minus,delta3minus,delta4minus]) = [40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,...
    -4,-4,-4,-4,-4,-4,-4,-4];
h=-w;

Inequality constraints;

A=zeros(22,51);
A(1,f89)=1;
A(2,f1011)=1;
A(3,f12fd)=1;
A(4,f13fd)=1;
A(5,f1314)=1;
A(6,f1215)=1;
A(7,f16fd)=1;
A(8,f17fd)=1;
A(9,f18fd)=1;
A(10,f19fd)=1;
A(11,f1920)=1;
A(12,f21fd)=1;
A(13,f2122)=1;
A(14,f24fd)=1;
A(15,f2425)=1;
A(16,f2627)=1;
A(17,f28fd)=1;
A(18,f2829)=1;
A(19,f30fd)=1;
A(20,f3031)=1;
A(21,f32fd)=1;
A(22,f3233)=1;

RHS for Inequality constraints;

b=zeros(22,1);
b(1)=9208333.33333333;
b(2)=4550000;
b(3)=2858333.33333333;
b(4)=3691666.66666667;
b(5)=3616666.66666667;
b(6)=5233333.33333333;
b(7)=3391666.66666667;
b(8)=5550000;
b(9)=2391666.66666667;
b(10)=42316666.6666667;
b(11)=408333.333333333;
b(12)=12808333.3333333;
b(13)=2391666.66666667;
b(14)=9466666.66666667;
b(15)=1175000;
b(16)=2733333.33333333;
b(17)=191091666.666667;
b(18)=2550000;
b(19)=1250000;
b(20)=2591666.66666667;
b(21)=1025000;
b(22)=5483333.33333333;

Equality constraints;

Aeq=zeros(22,51);
%Mass balance for node j
Aeq(1,[f25,f35,f56])=[1,1,-1];
Aeq(2,[f16,f56,f67])=[1,1,-1];
Aeq(3,[f47,f67,f78])=[1,1,-1];
Aeq(4,[f78,f89,f810])=[1,-1,-1];
Aeq(5,[f810,f1011,f1012])=[1,-1,-1];
Aeq(6,[f1012,f1213,f1215,f1216,f12fd])=[1,-1,-1,-1,-1];
Aeq(7,[f1213,f1314,f13fd])=[1,-1,-1];
Aeq(8,[f1216,f1617,f16fd])=[1,-1,-1];
Aeq(9,[f1617,f1718,f17fd])=[1,-1,-1];
Aeq(10,[f1718,f1819,f18fd])=[1,-1,-1];
Aeq(11,[f1819,f1920,f1921,f1923,f19fd])=[1,-1,-1,-1,-1];
Aeq(12,[f1921,f2122,f21fd])=[1,-1,-1];
Aeq(13,[f1923,f2324,f2326])=[1,-1,-1];
Aeq(14,[f2324,f2425,f24fd])=[1,-1,-1];
Aeq(15,[f2326,f2627,f2628])=[1,-1,-1];
Aeq(16,[f2628,f2829,f2830,f2832,f28fd])=[1,-1,-1,-1,-1];
Aeq(17,[f2830,f3031,f30fd])=[1,-1,-1];
Aeq(18,[f2832,f3233,f32fd])=[1,-1,-1];
%Modeling turn-up flexibility in the production fields.
Aeq(19,[f16,delta1minus,delta1plus])=[1,1,-1];
Aeq(20,[f25,delta2minus,delta2plus])=[1,1,-1];
Aeq(21,[f35,delta3minus,delta3plus])=[1,1,-1];
Aeq(22,[f47,delta4minus,delta4plus])=[1,1,-1];

RHS for equality constraints LB and UB;

beq=zeros(22,1);
beq(19)= 22583333.3333333;
beq(20)= 148325000;
beq(21)= 112516666.666667;
beq(22)= 107041666.666667;
lb = zeros(size(variables));
ub = Inf(size(variables));
ub(f16)=22583333.3333333;
ub(f25)=108325000;
ub(f35)=112516666.666667;
ub(f47)=107041666.666667;

Linprog;

startTime = tic;
[x, fval] = linprog(h,A,b,Aeq,beq,lb,ub);
time = toc(startTime);
for d = 1:N
fprintf('%12.2f \t%s\n',x(d),variables{d})
end

GA;

Objective=@Obj;
startTime = tic;
options=gaoptimset('UseParallel',true,'CreationFcn',@gacreationlinearfeasible);
[x,fval]=ga(Objective,51,A,b,Aeq,beq,lb,ub,[],options);
time_ga_sequential = toc(startTime);
fprintf('Serial GA optimization takes %g seconds.\n',time_ga_sequential);
for d = 1:N
fprintf('%12.2f \t%s\n',x(d),variables{d})
end

fmincon;

Objective=@Obj; 
startTime = tic;
  options = optimoptions(@fmincon,'UseParallel',true,'PlotFcns',@optimplotx);
  [x,fval] = fmincon(Objective,Initial,A,b,Aeq,beq,lb,ub,[],options);
time = toc(startTime);

for d = 1:N
fprintf('%12.2f \t%s\n',x(d),variables{d})
end

Here "Initial" is a solution from linprog! Only of I specify this starting point solver give me the feasible result otherwise it fails:(

[Anton Nester]

Linprog results:

16687556.82   f16
103639094.18   f25
100811633.32   f35
 94645049.01   f47
204450727.50   f56
221138284.32   f67
315783333.33   f78
306575000.00   f810
302025000.00   f1012
  7308333.33   f1213
286625000.00   f1216
283233333.33   f1617
277683333.33   f1718
275291666.67   f1819
 15200000.00   f1921
217366666.67   f1923
 10641666.67   f2324
206725000.00   f2326
203991666.67   f2628
  3841666.67   f2830
  6508333.33   f2832
  9208333.33   f89
  4550000.00   f1011
  2858333.33   f12fd
  3691666.67   f13fd
  3616666.67   f1314
  5233333.33   f1215
  3391666.67   f16fd
  5550000.00   f17fd
  2391666.67   f18fd
 42316666.67   f19fd
 12808333.33   f21fd
  2391666.67   f2122
   408333.33   f1920
  9466666.67   f24fd
  1175000.00   f2425
  2733333.33   f2627
191091666.67   f28fd
  2550000.00   f2829
  1250000.00   f30fd
  2591666.67   f3031
  1025000.00   f32fd
  5483333.33   f3233
  5895776.51   delta1minus
 44685905.82   delta2minus
 11705033.35   delta3minus
 12396617.66   delta4minus
        0.00   delta1plus
        0.00   delta2plus
        0.00   delta3plus
        0.00   delta4plus

Ga results;

    1773994.21   f16
   71627861.42   f25
   16426815.12   f35
   31567127.51   f47
   88054676.54   f56
   89828670.75   f67
  121395798.25   f78
  114368272.34   f810
  110635762.28   f1012
    5979573.70   f1213
   98573837.35   f1216
   97255480.85   f1617
   96001562.51   f1718
   95667570.47   f1819
    9612465.24   f1921
   72101219.13   f1923
    7458360.29   f2324
   64642858.84   f2326
   63227926.38   f2628
    2656258.28   f2830
    4843705.50   f2832
    7027525.91   f89
    3732510.07   f1011
    2238994.21   f12fd
    2976657.73   f13fd
    3002915.96   f1314
    3843357.02   f1215
    1318356.50   f16fd
    1253918.34   f17fd
     333992.03   f18fd
   13701265.78   f19fd
    8572368.50   f21fd
    1040096.74   f2122
     252620.32   f1920
    6906663.89   f24fd
     551696.40   f2425
    1414932.46   f2627
   53998857.42   f28fd
    1729105.17   f2829
     736098.42   f30fd
    1920159.87   f3031
     516179.92   f32fd
    4327525.58   f3233
   89598337.25   delta1minus
   76697138.72   delta2minus
  144336086.33   delta3minus
  130271329.54   delta4minus
   68788998.12   delta1plus
          0.13   delta2plus
   48246234.78   delta3plus
   54796790.38   delta4plus

Awful, isn't it? @gacreationlinearfeasible works like never before)))))

Fmincon gives the same solution as Linprog because I set up initial range.

[Alan Weiss]

You should take seriously the suggestion in the documentation and NOT give bounds as linear inequalities. You might also want to try linprog with various algorithms if you don't believe the result.

I would never use ga for this kind of problem, only fmincon once you include nonlinear constraints.

When I tried to run your linprog code I got the error

 Undefined function or variable 'f18fd'.

Alan Weiss

MATLAB mathematical toolbox documentation

[Anton Nester]

Thank you for your comments.I have changed the model according to your recommendations but the "ga" now throws the error:

Could not find a feasible initial point.
Serial GA optimization takes 0.242058 seconds.
Index exceeds matrix dimensions.

Error in MassBalance2 (line 128)
fprintf('%12.2f \t%s\n',x(d),variables{d})
If you have problems with code please find attached files with the Model.

[Anton Nester]

The main problem are periods! I need to calculate the flow for each pipeline through time horizon with changing demand and production with respect to quality, pressures and storages. And calculate the Initial point for each iteration for fmincon will be difficult task because I always need to compute linprog in parallel and implement it results to the fmincon. Thus I really need to configure one of those solvers to give me feasible solution for every iteration. This way I'm trying to solve it with "ga" also to compare the results and achieve more robust answers!!

[Alan Weiss]

Sorry, Anton, I still do not understand what you are trying to do, or why you are dissatisfied. In addition, I am surprised that you seem to continue trying to use ga, when it is almost certainly not an appropriate solver for your problem. If you have a smooth set of nonlinear constraints, fmincon is the solver of choice.

Check that your linear constraints are feasible by running linprog first. Once you get that to work, then you can try including nonlinear constraints using fmincon, and using the linprog solution as the initial point, as suggested here. In fact, that documentation link might help you solve your problem. Please check it out.

I am sorry if this does not help enough, but I really don't understand in what way you are dissatisfied.

Alan Weiss

MATLAB mathematical toolbox documentation

[Anton Nester]

Thank you! I see the point.