Why quadprog fvalue3 & fmincon fvalue1 are not same in the example below?

clc
clear
close all
price=[44.30133121	43.12406916	42.31422055	41.67599213	42.07729061	42.88582271	44.38989514	45.84411461	47.17838781	48.66472142	49.57772213	49.9537966	49.33055149	48.28029342	47.13492903	46.72706493	47.45610769	50.11226391	50.60314012	49.52144286	47.93957865	46.24335397	43.33797528	40.28594902	38.0842816	36.80017942	36.30394883	36.25929615	36.72432769	38.6619497	41.16986843	44.12309007	46.18842079	46.69566671	46.94333393	46.9464285	46.86775644	46.70612372	46.23813628	45.4851	47.01292074	51.29090431	53.48025645	52.85008289	52.1403893	50.98058129	48.48306443	45.56344424	43.8483262	43.21341571	42.97801815	43.29793883];
price=price';
inFlow=[  70	60	50	42	38	33	28	24	22	22	18	18	18	18	18	16	20	16	16	16	16	16	16	16	16	16	16	16	16	19	19	28	36	70	68	132	138	244	226	240	274	228	200	242	242	180	238	220	200	156	170	120 ];
inFlow=inFlow';
x0 = [inFlow; zeros(size(inFlow))];
stor0 = 180; 
k1 = 0.4422; 
k2 = 132.837; 
MW2kW = 1000; 
C2A = 1; 
C = [C2A k1 k2 MW2kW]; 
N = length(inFlow);
% DefineConstraints
LB = [9.072*ones(N,1); 
    zeros(N,1)];       
UB = [85.3*ones(N,1); 
    Inf(N,1)];       
ot = ones(N,1);
b = -5*ot;
A = spdiags([-ot -ot],[0 N],N,N*2);

A2 = spdiags([ot ot -ot -ot],[0 N -1 N-1],N,N*2);
b2 = 900*ot;

A2(1,:) = [];
b2(1,:) = [];

A = [A; A2; -A2];
b = [b; b2; b2];

c = stor0 + cumsum(inFlow); 
b = [b; 333.71-c; -63.51+c];
s = -sparse(tril(ones(N)));
s = [s s];
A = [A; s; -s];

Aeq = ones(1,2*N);
beq = sum(inFlow);

X = sym('x',[2*N,1]);   
F = sym('F',[N,1]);     
P = sym('P',[N,1]);     
s0 = sym('s0');         
c = sym({'c1','c2','c3','c4'}.','real');   

TotFlow = X(1:N)+X(N+1:end);
S = cell(N,1);
S{1} = s0 + c(1)*(F(1)-TotFlow(1));
for ii = 2:N
    S{ii} = S{ii-1} + c(1)*(F(ii)-TotFlow(ii));
end
k = c(2)*([s0; S(1:end-1)]+ S)/2+c(3);
MWh = k.*X(1:N)/c(4);
R = P.*MWh;
totR = -sum(R);
totR = subs(totR,[c;s0;P;F],[C';stor0;price;inFlow]);
matlabFunction(totR,'vars',{X},'file','objFcn');
options = optimset('MaxFunEvals',Inf,'MaxIter',1000,...
    'Algorithm','interior-point','Display','iter');
tic
[x1, fval1] = fmincon(@objFcn,x0,A,b,Aeq,beq,LB,UB,[],options);

Your initial point x0 is not between bounds lb and ub; FMINCON
shifted x0 to strictly satisfy the bounds.

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0     105   -5.317510e+04    1.965e+03    1.576e+01
    1     210   -5.325739e+04    1.945e+03    1.553e+01    3.859e+00
    2     315   -5.311530e+04    1.927e+03    1.295e+01    4.100e+00
    3     420   -5.188753e+04    1.824e+03    1.117e+01    2.721e+01
    4     525   -5.062621e+04    1.720e+03    1.055e+01    2.801e+01
    5     630   -4.907289e+04    1.591e+03    9.796e+00    3.493e+01
    6     735   -4.614195e+04    1.345e+03    9.526e+00    6.703e+01
    7     840   -4.352817e+04    1.122e+03    8.785e+00    6.112e+01
    8     945   -4.323210e+04    1.097e+03    8.785e+00    6.820e+00
    9    1050   -4.170045e+04    9.668e+02    8.785e+00    3.630e+01
   10    1155   -3.766727e+04    6.258e+02    8.785e+00    9.469e+01
   11    1260   -3.633916e+04    5.246e+02    8.785e+00    2.706e+01
   12    1365   -3.313963e+04    2.779e+02    3.875e+00    6.683e+01
   13    1470   -3.202658e+04    1.886e+02    3.875e+00    2.580e+01
   14    1575   -3.181577e+04    1.698e+02    3.875e+00    7.178e+00
   15    1680   -3.169923e+04    1.588e+02    3.875e+00    5.054e+00
   16    1785   -3.160061e+04    1.488e+02    3.875e+00    5.994e+00
   17    1890   -3.156994e+04    1.405e+02    3.875e+00    1.249e+01
   18    1995   -3.153315e+04    1.376e+02    3.875e+00    8.783e-01
   19    2100   -3.126567e+04    1.149e+02    4.068e+00    8.899e+00
   20    2205   -3.105143e+04    9.704e+01    3.175e+00    7.159e+00
   21    2310   -3.093686e+04    8.575e+01    3.030e+00    8.834e+00
   22    2415   -3.088176e+04    8.088e+01    3.039e+00    2.584e+00
   23    2520   -3.071320e+04    6.663e+01    3.056e+00    6.251e+00
   24    2625   -3.058859e+04    5.576e+01    3.075e+00    6.054e+00
   25    2730   -3.053859e+04    5.081e+01    3.095e+00    4.842e+00
   26    2835   -3.048442e+04    4.582e+01    2.130e+00    3.994e+00
   27    2940   -3.007487e+04    1.409e+01    2.148e+00    9.547e+00
   28    3045   -3.005304e+04    1.220e+01    2.151e+00    1.440e+00
   29    3150   -3.001689e+04    9.048e+00    2.001e+00    3.123e+00
   30    3255   -2.999671e+04    7.347e+00    2.001e+00    1.292e+00

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   31    3360   -2.998759e+04    5.995e+00    2.001e+00    4.888e+00
   32    3465   -2.998729e+04    5.522e+00    2.001e+00    3.558e+00
   33    3570   -2.998258e+04    4.695e+00    1.983e+00    3.889e+00
   34    3675   -2.997754e+04    3.553e+00    1.943e+00    7.068e+00
   35    3780   -2.997669e+04    3.405e+00    1.937e+00    9.695e-01
   36    3885   -2.997976e+04    3.288e+00    2.067e+00    4.067e+00
   37    3990   -2.998450e+04    3.212e+00    2.067e+00    5.174e+00
   38    4095   -2.998058e+04    2.909e+00    1.793e+00    1.323e-01
   39    4200   -2.997060e+04    2.060e+00    1.646e+00    1.130e+00
   40    4305   -2.994842e+04    0.000e+00    4.340e-01    6.258e+00
   41    4410   -2.994965e+04    0.000e+00    4.175e-01    1.877e+00
   42    4515   -2.994983e+04    0.000e+00    4.149e-01    2.876e-01
   43    4620   -2.995022e+04    0.000e+00    3.179e-01    7.384e-01
   44    4725   -2.995205e+04    9.095e-13    3.331e-01    4.167e+00
   45    4830   -2.995526e+04    0.000e+00    3.161e-01    8.108e+00
   46    4935   -2.995696e+04    0.000e+00    3.064e-01    5.050e+00
   47    5041   -2.995697e+04    0.000e+00    2.735e-01    2.474e-02
   48    5146   -2.995930e+04    0.000e+00    3.033e-01    8.306e+00
   49    5251   -2.996097e+04    9.095e-13    3.999e-01    8.181e+00
   50    5356   -2.996234e+04    0.000e+00    4.397e-01    1.213e+01
   51    5461   -2.996236e+04    0.000e+00    8.303e-02    1.476e-01
   52    5566   -2.996265e+04    0.000e+00    6.168e-02    1.042e+00
   53    5671   -2.996276e+04    0.000e+00    9.242e-02    2.505e+00
   54    5776   -2.996277e+04    9.095e-13    4.885e-02    5.850e-01
   55    5881   -2.996278e+04    0.000e+00    4.551e-02    2.602e-01
   56    5986   -2.996281e+04    0.000e+00    4.597e-02    6.206e-01
   57    6091   -2.996297e+04    9.095e-13    4.561e-02    2.979e+00
   58    6196   -2.996369e+04    0.000e+00    5.375e-02    1.461e+01
   59    6301   -2.996457e+04    0.000e+00    1.099e-01    1.891e+01
   60    6407   -2.996457e+04    0.000e+00    1.099e-01    4.754e-02

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   61    6512   -2.996458e+04    0.000e+00    4.863e-02    3.276e-01
   62    6617   -2.996459e+04    0.000e+00    5.321e-02    1.302e+00
   63    6722   -2.996460e+04    0.000e+00    5.923e-02    5.998e-01
   64    6827   -2.996466e+04    0.000e+00    8.090e-02    3.031e+00
   65    6932   -2.996476e+04    9.095e-13    1.000e-01    4.534e+00
   66    7037   -2.996478e+04    0.000e+00    8.306e-02    1.477e+00
   67    7142   -2.996485e+04    9.095e-13    1.952e-02    2.486e+00
   68    7247   -2.996492e+04    0.000e+00    8.790e-03    8.749e-01
   69    7352   -2.996492e+04    9.095e-13    7.346e-04    2.524e-02
   70    7457   -2.996493e+04    0.000e+00    8.262e-04    1.961e-02
   71    7562   -2.996493e+04    0.000e+00    6.363e-04    3.319e-03
   72    7667   -2.996493e+04    0.000e+00    6.349e-04    9.708e-04
   73    7785   -2.996493e+04    0.000e+00    1.190e-05    9.091e-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.

fval1

fval1 = -2.9965e+04

toc

Elapsed time is 1.872455 seconds.

Rsub = subs(R,[c;s0;P;F],[C';stor0;price;inFlow]);
Rtot = -sum(Rsub);
fsym = gradient(Rtot,X);
f = double(subs(fsym,X,zeros(size(X))));
H = 0.5.*double(hessian(Rtot,X));
qpoptions = optimset('Algorithm','interior-point-convex','Disp','iter');

tic
[x3,fval3] = quadprog(H,f,A,b,Aeq,beq,LB,UB,[],qpoptions);

Solved 0 variables, 0 equality, and 55 inequality constraints during the presolve.

 Iter            Fval  Primal Infeas    Dual Infeas  Complementarity  
    0   -1.124215e+05   8.371521e+04   1.889709e+03     9.448545e+02  
    1   -1.112268e+05   8.361616e+04   1.888309e+03     9.440566e+02  
    2   -1.057057e+05   7.905741e+04   1.787875e+03     8.973915e+02  
    3   -8.129912e+04   5.760112e+04   1.319815e+03     6.265520e+02  
    4   -7.868598e+04   9.299590e+03   2.064001e+02     9.231524e+01  
    5   -5.210382e+04   4.432996e+00   1.032001e-01     3.238543e+01  
    6   -6.699990e+04   1.218708e-01   2.837150e-03     1.147093e+01  
    7   -7.193226e+04   4.069387e-02   9.473523e-04     6.299708e+00  
    8   -7.432700e+04   2.143778e-03   4.990709e-05     3.509479e+00  
    9   -7.458784e+04   3.248514e-04   7.562533e-06     1.253393e+00  
   10   -7.465602e+04   2.979443e-05   6.936132e-07     1.435589e-01  
   11   -7.466368e+04   6.509312e-07   1.515396e-08     2.237764e-02  
   12   -7.466409e+04   2.144702e-09   5.420900e-11     1.406085e-03  
   13   -7.466409e+04   4.365575e-11   2.129399e-13     4.247797e-07  
   14   -7.466409e+04   2.776801e-11   2.842171e-14     1.852086e-17  

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.

fval3

fval3 = -7.4664e+04

toc

Elapsed time is 0.170306 seconds.

norm(x1-x3)

ans = 153.3015

Answer 1

[Torsten]

Should be

H = double(hessian(Rtot,X));

instead of

H = 0.5*double(hessian(Rtot,X));

since

Rtot(x) = Rtot(0) + f.'*x + 0.5*x.'*H*x

So you solve two different problems with two different solutions.

[Chandra Upadhyaya]

The method suggested by is not acceptable. I have taken reference to write the code from "Optimization in MATLAB: An Introduction to Quadratic Programming" by Seth DeLand (2012). He is the marketing manager from MATLAB. I am supposed to get both fvalue1 & fvalue3 same. Also x1 and x3 should be same.

[Torsten]

I am supposed to get both fvalue1 & fvalue3 same. Also x1 and x3 should be same.

No.

The Taylor series development of Rtot is as written above. Thus in order to minimize Rtot, you will have to supply H = double(hessian(Rtot,X)) to "quadprog", not H = 0.5*double(hessian(Rtot,X)).

Make the test:

If you use "fmincon" to solve the problem for which you use "quadprog" at the moment, you will get the same result as quadprog gives.