How to avoid a negative solution with fmincon, including an external equation?

Hello community,

I want to create an optimization script with fmincon, including an equation (f2), which is not solved by fmincon. This equation will be solved with a function (see bottom of the script in the function 'kostfunktion'), but it contains the parameters x(1) and x(2) from the optimization problem. The problem is, that f2 and therefore fval are always negative. How can I set these parameters to a minimum of 0 or at least let them be positive?

The involving case is this one:

• f1 = 10*x1 + 2*x2 - 30 f2 = - 2*x1 - 10*x2 + 30

Bounds: 1<=x1<=5 1<=x2<=5 0<=f1<=inf 0<=f2<=inf

The optimization goal is to get the minimum of: obj = f1 + f2

f1 will be solved by fmincon. f2 will be solved in the function 'kostfunktion'.

clc, clear, clear global;
global f2   %f2 appears in the workspace, just to be visible
x0 = [2; 2; 1]; % [x1, x2, f1]
A = [];
b = [];
Aeq = [-10 -2 1];
beq = [-30];
lb = [1; 1; 0];
ub = [5; 5; inf];
options = optimoptions('fmincon', 'Display', 'iter');
[x fval] = fmincon(@kostfunktion,x0,A,b,Aeq,beq,lb,ub,[],options);

function obj = kostfunktion(x)
global f2
f2 = -2*x(1)-10*x(2)+30;
    obj = x(3) + f2;    % obj = f1 (solved with fmincon) + f2 (solved with this function)
end

The optimal solution should be obj = 0 with x1 = 2.5 and x2 = 2.5, but in this script, the solution is obj = -24, because I can't set bounds for f2. How can I set the bounds for f2 or fval?

Thanks for your support.

Answer 1

[Bruno Luong]

The condition

f2 = - 2*x1 - 10*x2 + 30 >= 0

is equivalen to

2*x1 + 10*x2 + 0*f1 <= 30.

So set

A = [2 10 0]
b = 30

If you don't want to use A/b arguments (since f2 can be changed to more generic non linear function), just use non linear constraint nonlcon argument and program c output that handles that constraint.

[Drilon]

Thanks Bruno,

I will test it.

[Matt J]

Why is it okay to change A,b but not Aeq,beq ...?

[Drilon]

I think he did not wrote that I have to change A,b. He is suggesting a new program c which handles the constraint. I will try it out and write about the process.

[Bruno Luong]

Use nonlcon (less performance if you apply to a linear constraint than setting A/b)

x0 = [2; 2; 1]; % [x1, x2, f1]
A = [];
b = [];
Aeq = [-10 -2 1];
beq = [-30];
lb = [1; 1; 0];
ub = [5; 5; inf];
options = optimoptions('fmincon', 'Display', 'iter');
[x fval] = fmincon(@kostfunktion,x0,A,b,Aeq,beq,lb,ub,@mycon,options);

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       4    7.000000e+00    7.000e+00    3.994e+00
    1       8    6.336764e-01    2.737e+00    1.562e+00    7.264e-01
    2      12    1.155960e-01    1.174e-02    4.394e-01    6.429e-01
    3      16    2.008181e-02    1.676e-03    2.006e-02    1.002e-02
    4      20    1.993555e-03    3.553e-15    1.003e-03    2.173e-03
    5      24    2.051120e-05    0.000e+00    1.907e-04    1.007e-03
    6      28    1.999992e-05    0.000e+00    1.000e-05    5.092e-07

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

x

x = 3×1

2.5000 2.5000 0.0000

function obj = kostfunktion(x)
global f2
f2 = -2*x(1)-10*x(2)+30;
    obj = x(3) + f2;    % obj = f1 (solved with fmincon) + f2 (solved with this function)
end

function [c,ceq] = mycon(x)
f2 = -2*x(1)-10*x(2)+30;
c = -f2; % (c <= 0) <=> (f2 >= 0)
ceq = [];
end

[Bruno Luong]

Use A/b

x0 = [2; 2; 1]; % [x1, x2, f1]
A = [2 10 0]

A = 1×3

2 10 0

b = 30;
Aeq = [-10 -2 1];
beq = [-30];
lb = [1; 1; 0];
ub = [5; 5; inf];
options = optimoptions('fmincon', 'Display', 'iter');
[x fval] = fmincon(@kostfunktion,x0,A,b,Aeq,beq,lb,ub,[],options);

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       4    7.000000e+00    7.000e+00    3.994e+00
    1       8    6.336764e-01    2.737e+00    1.562e+00    7.264e-01
    2      12    1.155960e-01    1.174e-02    4.394e-01    6.429e-01
    3      16    2.008181e-02    1.676e-03    2.006e-02    1.002e-02
    4      20    1.993555e-03    0.000e+00    1.003e-03    2.173e-03
    5      24    2.051125e-05    0.000e+00    1.908e-04    1.007e-03
    6      28    1.999992e-05    3.553e-15    1.000e-05    5.092e-07

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

x

x = 3×1

2.5000 2.5000 0.0000

function obj = kostfunktion(x)
global f2
f2 = -2*x(1)-10*x(2)+30;
    obj = x(3) + f2;    % obj = f1 (solved with fmincon) + f2 (solved with this function)
end

[Drilon]

Thanks Bruno, this solution works too.

I will try to understand your solution too and test it for another case.

Answer 2

[Matt J]

Since your problem is linear, let's just use linprog,

%Reformulate with x=[x1,x2,f1,f2]

A = [];
b = [];
Aeq = [-10 -2  1  0; 
        2  10  0  1 ];
beq = [-30;+30];
lb = [1; 1; 0; 0];
ub = [5; 5; inf; inf];

obj=[0 0 1 1];
x=linprog(obj,A,b,Aeq,beq,lb,ub)

Optimal solution found.

x = 4×1

2.5000 2.5000 0 0

[Drilon]

Hello Matt,

thanks for the answer.

I need to solve f2 without fmincon, since this function will be looped later (it's a parameter). The equation 'Aeq*x=beq' in linprog can not be changed afterward, therefore I need to solve the function f2 without fmincon.

[Matt J]

I did the whole thing for you without fmincon!

And of course you can modify Aeq,beq in a loop.

[Drilon]

Yeah, but I don't need linprog. This script is just a testing script for the objective script, which includes much more.

And no, you can't modify Aeq, beq in a loop, since the changes won't apply for the equation 'Aeq*x=beq'.

The goal is, to include an equation like f2, without entering it in the optimization parameters.

[Matt J]

And no, you can't modify Aeq, beq in a loop

You have to tell us what we are free to modify... Only the objective function code?

[Drilon]

Dude, I wrote that you can't modify it at all, regarding to Mathworks. If you can do it, you could show me how. Because whenever I do it, the optimization solver still uses the old numbers.

The easiest solution would be the modification of beq, but it doesn't work, Mathworks told me. If you know how, we would be happy about the information.

[Matt J]

Something has to be changed in order to produce different results....

[Drilon]

Yes but regarding to the Mathworks Service, I can't change the constraints. I tried it and my code was not able to change it, because the optimization continued with the old constraints.

[Matt J]

I don't understand what the "MathWorks Service" refers to here. You said you have many optimization problems to solve and you are looping over them, each time calling fmincon to get a solution. Why can't you feed a different A,b,Aeq,beq to fmincon each time you call it?

[Drilon]

Missunderstanding. I have only one optimization problem with different constraints. The constraints will change for many cases and all these cases need to fit into one optimization. There will be million cases, therefore I do not want to do it without parameters, I need them to create call all the cases.

Answer 3

[Matt J]

You could consider something like this:

x0 = [2; 2; 1]; % [x1, x2, f1]
A = [];
b = [];
Aeq = [-10 -2 1];
beq = [-30];
lb = [1; 1; 0];
ub = [5; 5; inf];
options = optimoptions('fmincon', 'Display', 'iter','Algorithm', 'sqp','MaxFunctionEvaluations',inf,...
                       'OptimalityTolerance',1e-12,'FunctionTolerance',1e-12,'StepTolerance',1e-12);
[x fval] = fmincon(@kostfunktion,x0,A,b,Aeq,beq,lb,ub,[],options)

 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0           4    7.000000e+00     7.000e+00     1.000e+00     0.000e+00     1.000e+01  
    1          13    1.816455e+00     5.824e+00     1.681e-01     5.219e-01     8.824e+00  
    2          26    1.018196e+00     5.589e+00     4.035e-02     1.066e-01     7.560e+01  
    3          30    8.273146e-01     3.553e-15     1.000e+00     9.951e-01     1.012e+01  
    4          34    7.326886e-01     3.553e-15     1.000e+00     1.005e-02     1.001e+01  
    5          38    2.595588e-01     0.000e+00     1.000e+00     5.026e-02     1.001e+01  
    6          48    1.874546e-02     0.000e+00     1.176e-01     2.957e-02     2.062e+01  
    7          57    1.274881e-02     0.000e+00     1.681e-01     3.347e-03     1.001e+01  
    8          61    2.463475e-03     3.553e-15     1.000e+00     1.093e-03     1.000e+01  
    9          70    1.949447e-03     0.000e+00     1.681e-01     9.182e-04     6.375e+00  
   10          74    3.999204e-04     3.553e-15     1.000e+00     3.575e-04     2.855e+00  
   11          78    3.656302e-07     0.000e+00     1.000e+00     2.899e-04     8.573e-02  
   12          82    5.079048e-11     0.000e+00     1.000e+00     8.976e-06     1.200e-03  
   13          86    5.079048e-11     0.000e+00     1.000e+00     1.939e-21     1.108e-03  

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

x = 3×1

2.5000 2.5000 0.0000

fval = 5.0790e-11

function obj = kostfunktion(x)
global f2


    f2 =-2*x(1)-10*x(2)+30;

   c=10;
   if f2<0, f2=exp(-c*f2)+c*f2-1; end %apply a penalty

    obj = x(3) + f2;    % obj = f1 (solved with fmincon) + f2 (solved with this function)
end

[Drilon]

Thanks Matt, it works with this code.

I will try to figure out how the 'if' expression (the penalty) works and how it affects my solving process.

Could you explain it shortly?

[Matt J]

The penalty pushes the objective steeply up toward infinitey as soon as f2 tries to become negative. Therefore, fmincon is discouraged from looking for a solution there.

[Drilon]

That means, fmincon tries to find another solution.

I can see that the code works, but I would like to ask: can there be any negative consequences for the further steps or the final solution, since it is pushing f2 toward infinitey?

[Matt J]

Well in general this is an approximate solution, one which gets better and better as the steepness parameter c is increased.

I would much prefer to handle this through the constraints, but from your earlier comments it is not clear what freedom you have to change those. You mentioned that these optimizations are occuring in a loop during which constraint data like Aeq, beq is forbidden to vary.

[Drilon]

Dear Matt, please have a look at the following code. I just found it. It contains the described problem, where you can't change the value of beq, because fmincon is not able to change the value after defining it.

However, your solution on the top works well. I will try it for another case and close this topic if everything works fine.

clc, clear, clear global;
global Variabel2
if exist(Variabel2)==0 %falls Variabel2 nicht existiert, wird diese erstellt
        Variabel2=20; %Variabel2 is declared as 20
end
%% Optimierung durch 'fmincon'
x0 = [1; 1; 0; 0];
A = [];
b = [];
Aeq = [-10 -2 1 0; 2 10 0 1];
global beq
beq = [-30; Variabel2]; %beq übernimmt Variabel2 als Wert
lb = [1; 1; 0; 0];
ub = [5; 5; inf; inf];
options = optimoptions('fmincon', 'Display', 'iter');

[x fval] = fmincon(@kostfunktion,x0,A,b,Aeq,beq,lb,ub,[],options);

function obj = kostfunktion(x)
    global Variabel2
    Variabel2=30;   
    global beq
    beq(2,1)=Variabel2      %Variabel2 is now placed in beq with the new value 30
    % this value can be seen in the workspace, but it has no effect at all, because
    % the optimization continues with the first definition of beq which was
    % beq = [-30;20]. You can see that on the results. x(1) ist 2.7 and x(2) is 1.45.
    % If you change Variabel2 under the 'exist' command to 30, you will get
    % the optimal solution with 30.
    obj = x(3) + x(4);
end

[Bruno Luong]

"It contains the described problem, where you can't change the value of beq, because fmincon is not able to change the value after defining it. "

This code does something definitivelty not allowed and ineffective (global variable is NOT in the scope of FMINCON). The author of this code does not know program correctly MATLAB.

But Matt never did something like this.

I think YOU just miss-read his posts.

[Matt J]

@Drilon

Yes, but why would you be attempting to make beq(2) evolve throughout the optimization?

If you have a constraint of the general form Aeq*x=h(x) where h(x) is a non-constant, non-linear function of x, then this is a non-linear constraint, and should be handled using the nonlcon input.

[Drilon]

Bruno, it was not his code. I showed him an old code where the problem appeared, because he wanted to understand the problem about beq. Because of this problem, I was looking for another solution, for example like the suggestion of Matt or yours.

However, thanks guys. Both solutions are great, but I think that nonlcon fits the most for my problem.

Thanks for the help. I appreciate it :)