fmincon get the wrong answer

i have a optimization problem that fmincon get the wrong answer, there are a demand power(2) that can be generate or buy with different price, the codes are as follow:

a = 0.005;
b = 6;
c = 100;

% Define the cost function for generating energy (quadratic cost function) for the current microgrid
Cg_i = @(Ec, a, b, c) a * Ec^2 + b * Ec + c;

% Given prices
prices = [33, 26];
mprice = 27; % Given value of mprice

% Define the cost function for buying energy
Cb_i = @(Ec, price) price * Ec;

% Given total energy constraint
Ec_total = 2; % Given value of Ec_total

% Define the objective function
objective = @(X) Cg_i(X(1), a, b, c) + Cb_i(X(2), prices(1)) + Cb_i(X(3), prices(2)) + Cb_i(X(4), mprice);

% Initial guess for Ec_g, Ec1, Ec2
initial_guess = [Ec_total / 3, Ec_total / 3, Ec_total / 3, Ec_total / 3];
% initial_guess = [0, 0, Ec_total, 0];
% Linear equality constraint for the sum of Ec_g, Ec1, Ec2
Aeq = [1, 1, 1, 1];
beq = Ec_total;

% Options for fmincon
options = optimoptions('fmincon', 'Display', 'iter');
% Define solver options

% Call fmincon with specified options
X = fmincon(objective, initial_guess, [], [], Aeq, beq, [0, 0, 0,0], [], [], options)

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       5    1.613356e+02    6.667e-01    1.700e+01
    1      10    1.244878e+02    5.104e-01    1.585e+01    1.761e+00
    2      15    1.175907e+02    4.127e-01    6.932e-01    2.443e-01
    3      20    1.122426e+02    1.234e-03    1.529e-01    3.250e-01
    4      25    1.120220e+02    2.220e-16    3.175e-03    1.134e-02
    5      30    1.120200e+02    0.000e+00    5.612e-06    1.112e-04

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 = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

the matlab answer is X=[2 0 0 0]; that means generat 2 and no buy, that obviously is wrong because cost of buy with price 26 is 26*2=52 while generating cost is 112!

Answer 1

[Torsten]

Cost of buy with price 26 is 2*26 + c = 152 in your setting. Note that although you don't generate energy, Cg_i(0,a,b,c) = c - thus you always add c to the cost of buying energy which obviously is not correct. Or do you want to treat c as fixed costs that always have to be paid - independent of whether energy is generated or not ?

[encieh baqeri]

thanks, yes i didnt see that! the constant c should be add when the generating select not always! do you have any suggestion for this?

[Torsten]

One way is to use "ga": it can handle discontinuous objective functions:

rng("default")

a = 0.005;
b = 6;
c = 100;

% Define the cost function for generating energy (quadratic cost function) for the current microgrid
Cg_i = @(Ec, a, b, c) (a * Ec^2 + b * Ec + c)*(Ec>0);

% Given prices
prices = [33, 26];
mprice = 27; % Given value of mprice

% Define the cost function for buying energy
Cb_i = @(Ec, price) price * Ec;

% Given total energy constraint
Ec_total = 2; % Given value of Ec_total

% Define the objective function
objective = @(X) Cg_i(X(1), a, b, c) + Cb_i(X(2), prices(1)) + Cb_i(X(3), prices(2)) + Cb_i(X(4), mprice);

% Initial guess for Ec_g, Ec1, Ec2
initial_guess = [Ec_total / 3, Ec_total / 3, Ec_total / 3, Ec_total / 3];
% initial_guess = [0, 0, Ec_total, 0];
% Linear equality constraint for the sum of Ec_g, Ec1, Ec2
Aeq = [1, 1, 1, 1];
beq = Ec_total;

X = ga(objective,4,[],[],Aeq,beq,zeros(4,1),inf(4,1))

ga stopped because the average change in the fitness value is less than options.FunctionTolerance.

X = 1x4

0 0.0000 2.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

objective(X)

ans = 52.0000

[Torsten]

Maybe somebody sees the problem - I don't know why "fmincon" has to shift the initial guess to [1 1 2 1] . In my opinion, [0 0 2 0] as given should satisfy lower and upper bounds as well as the linear equality constraint.

% Given total energy constraint
Ec_total = 2; % Given value of Ec_total
% Initial guess for Ec_g, Ec1, Ec2
initial_guess = [0, 0, Ec_total, 0];
% Linear equality constraint for the sum of Ec_g, Ec1, Ec2
Aeq = [1, 1, 1, 1];
beq = Ec_total;
% Options for fmincon
options = optimoptions('fmincon', 'Display', 'iter');
% Define solver options
% Call fmincon with specified options
[X,fval,flag] = fmincon(@objective, initial_guess, [], [], Aeq, beq, zeros(1,4), inf(1,4), [], options)

Initial point X0 is not between bounds LB and UB; 
FMINCON shifted X0 to strictly satisfy the bounds.

X = 1x4

0.9900 0.9900 2.0000 0.9900

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

0.9900 0.9900 2.0000 0.9900

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

0.9900 0.9900 2.0000 0.9900

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

0.9900 0.9900 2.0000 0.9900

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

0.9900 0.9900 2.0000 0.9900

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       5    2.173449e+02    2.970e+00    1.739e+01

X = 1x4

3.8590 0.2305 0.0100 0.0220

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

3.8590 0.2305 0.0100 0.0220

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

3.8590 0.2305 0.0100 0.0220

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

3.8590 0.2305 0.0100 0.0220

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

3.8590 0.2305 0.0100 0.0220

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

    1      10    1.316905e+02    2.122e+00    1.666e+01    3.702e+00

X = 1x4

3.8050 0.1582 0.0001 0.0112

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

3.8050 0.1582 0.0001 0.0112

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

3.8050 0.1582 0.0001 0.0112

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

3.8050 0.1582 0.0001 0.0112

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

3.8050 0.1582 0.0001 0.0112

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

    2      15    1.284243e+02    1.974e+00    1.087e+00    9.153e-02

X = 1x4

2.0519 0.0008 0.0048 0.0049

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0519 0.0008 0.0048 0.0049

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0519 0.0008 0.0048 0.0049

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0519 0.0008 0.0048 0.0049

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0519 0.0008 0.0048 0.0049

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

    3      20    1.126136e+02    6.233e-02    5.803e-01    1.760e+00

X = 1x4

1.9972 0.0007 0.0010 0.0010

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

1.9972 0.0007 0.0010 0.0010

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

1.9972 0.0007 0.0010 0.0010

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

1.9972 0.0007 0.0010 0.0010

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

1.9972 0.0007 0.0010 0.0010

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

    4      25    1.120824e+02    1.110e-15    1.078e-02    5.499e-02

X = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

    5      30    1.120206e+02    2.220e-16    6.481e-04    3.221e-03

X = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

    6      35    1.120200e+02    0.000e+00    9.988e-07    3.269e-05

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 = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fval = 112.0200

flag = 1

objective(X)

X = 1x4

2.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = 112.0200

function value = objective(X)
  X
  a = 0.005;
  b = 6;
  c = 100;
  % Define the cost function for generating energy (quadratic cost function) for the current microgrid
  Cg_i = @(Ec, a, b, c) a * Ec^2 + b * Ec + c;
  % Given prices
  prices = [33, 26];
  mprice = 27; % Given value of mprice
  % Define the cost function for buying energy
  Cb_i = @(Ec, price) price * Ec;
  % Define the objective function
  value = Cg_i(X(1), a, b, c) + Cb_i(X(2), prices(1)) + Cb_i(X(3), prices(2)) + Cb_i(X(4), mprice);
end