Surface plot doesn't correspond to optimization solution

My problem is that the surface plot doesn't correspond to values obtained from the optimizaiton problem.

I want to plot an input mix over different calibrations for two parameters. Hence, I solve the optimization problem in two different loops, see below, and store the results in several matrices. The problem is now, when I select a point randomly in the surface plot, and then calibrate the model for these values, I do not obtain the same solution, Xop(3). What is the reason for this? Any suggestion is highly appreciated, thank you.

%% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP;
alphaK = 0.6;
alphaL = 0.4;
bbeta = 0.8;
delta = 0.4;
%phi = 0.89;
rho = 0.19;
%v = 0.0;
w = 1;
mu = 1;                                % Mean
sdev = 0.1;                            % Standard deviation
rng(0,'twister')
p = mu + sdev.*randn(10000,1);
cevP = 1;

PHI = [0:0.01:1];                      % Set \phi and \nu values that should vary.
NU = [0:0.0001:0.01];

K = zeros(size(PHI,2), size(NU,2));    % Empty matrices for loop
B = zeros(size(PHI,2), size(NU,2));
Fp = zeros(size(PHI,2), size(NU,2));
MAX = zeros(size(PHI,2), size(NU,2));

%% Solve Optimization problem in Loop over different values
for i = 1:size(PHI, 2)
   
    phi = PHI(i);
    
    for j = 1:size(NU,2)

        v = NU(j);

x0 = [2 2 0.5 1];                      % Start values for optimization; K, B, CDF, p_bar
lb = [0 0 0 0];                        % Lower bound
ub = [Inf Inf 1 Inf];                  % Upper bound
nonlcon = @constraint;
Objfcn = @optmprob;
options = optimoptions(@fmincon,'Algorithm','interior-point','Display','iter');
[Xop, Fop] = fmincon(Objfcn,x0,[],[],[],[],lb,ub,nonlcon,options)
K(i,j) = Xop(1);
B(i,j) = Xop(2);
Fp(i,j) = Xop(3);
MAX(i,j) = -Fop;
end
end
L = (B + NU)./ w;                      % Define L
Imix = L ./ (L+K);                     % Define input

figure()
[X, Y] = meshgrid(PHI, NU);
surf(X,Y,Fp)
xlabel('\phi')
ylabel('\nu')
zlabel('F')


function profit = optmprob(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP;
Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
cevP = quadgk(@(p) (p.*normpdf(p,mu,sdev))./(1-x(3)), x(4), Inf);       % Conditional expected value in integral form

profit = -(x(3).*phi.*(v+x(2)+delta.*x(1)) + cevP.*Yprod - (x(1)+x(2)));
end


function [c,ceq] = constraint(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev p_bar;
Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);

ceq(1) = x(3) - normcdf(x(4),mu,sdev);
ceq(2) = x(4) - ((phi.*(v+x(2)+delta.*x(1)) + (1+(x(3)./(1-x(3)))).*(x(1)+x(2)))./Yprod);
c = [];
end

Answer 1

[Torsten]

The problem is now, when I select a point randomly in the surface plot, and then calibrate the model for these values, I do not obtain the same solution, Xop(3).

I don't know what you mean.

If you mean the error message that the dimensions for the surface plot don't match:

You must use

surf(X,Y,Fp.')

instead of

surf(X,Y,Fp)

[TerribleStudent]

Yes that is exactly it, thank you. So my question would be, how did you identify the issue? Should I know this from simple linear algebra? By just running my code above, doesn't give any error messages etc. Sorry for this rather novice question. Thank you.

[Torsten]

So my question would be, how did you identify the issue? Should I know this from simple linear algebra?

No, you just have to read the documentation of surf:

surf(X,Y,Z)

X: x-coordinates, specified as a matrix the same size as Z, or as a vector with length n, where [m,n] = size(Z).

Y: y-coordinates, specified as a matrix the same size as Z, or as a vector with length m, where [m,n] = size(Z).

Z: z-coordinates, specified as a matrix.

Thus in your code, you can either change the roles of X and Y or you can transpose Z in the call to surf.

But I must admit that the inputs for surface plots confuse me, too, every time I use them. Same for the multidimensional interpolation routines.

[TerribleStudent]

I'm still a bit confused, but thank you.

[Torsten]

The number of rows of Z must be equal to the length of the Y-vector (not the X-vector !), the number of columns of Z must be equal to the length of the X-vector (not the Y-vector !).

[TerribleStudent]

Thanks. So in principle I could have "filled" the Z matrix Z(j,i) instead of Z(i,j).

What about the two spikes in the plot? Is it common for optimization problems?

[Torsten]

I cannot reproduce spikes because computation time with MATLAB online is too short.

%% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP;
alphaK = 0.6;
alphaL = 0.4;
bbeta = 0.8;
delta = 0.4;
%phi = 0.89;
rho = 0.19;
%v = 0.0;
w = 1;
mu = 1;                                % Mean
sdev = 0.1;                            % Standard deviation
rng(0,'twister')
p = mu + sdev.*randn(10000,1);
cevP = 1;

PHI = [0:0.1:1];                      % Set \phi and \nu values that should vary.
NU = [0:0.001:0.01];

K = zeros(size(PHI,2), size(NU,2));    % Empty matrices for loop
B = zeros(size(PHI,2), size(NU,2));
Fp = zeros(size(PHI,2), size(NU,2));
MAX = zeros(size(PHI,2), size(NU,2));

%% Solve Optimization problem in Loop over different values
for i = 1:size(PHI, 2)
   
    phi = PHI(i);
    
    for j = 1:size(NU,2)

        v = NU(j);

x0 = [2 2 0.5 1];                      % Start values for optimization; K, B, CDF, p_bar
lb = [0 0 0 0];                        % Lower bound
ub = [Inf Inf 1 Inf];                  % Upper bound
nonlcon = @constraint;
Objfcn = @optmprob;
options = optimoptions(@fmincon,'Algorithm','interior-point','Display','none');
[Xop, Fop] = fmincon(Objfcn,x0,[],[],[],[],lb,ub,nonlcon,options);
K(i,j) = Xop(1);
B(i,j) = Xop(2);
Fp(i,j) = Xop(3);
MAX(i,j) = -Fop;
end
end
L = (B + NU)./ w;                      % Define L
Imix = L ./ (L+K);                     % Define input

figure()
[X, Y] = meshgrid(PHI, NU);
surf(X,Y,Fp.')
xlabel('\phi')
ylabel('\nu')
zlabel('F')

function profit = optmprob(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP;
Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
cevP = quadgk(@(p) (p.*normpdf(p,mu,sdev))./(1-x(3)), x(4), Inf);       % Conditional expected value in integral form
profit = -(x(3).*phi.*(v+x(2)+delta.*x(1)) + cevP.*Yprod - (x(1)+x(2)));
end
function [c,ceq] = constraint(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev p_bar;
Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
ceq(1) = x(3) - normcdf(x(4),mu,sdev);
ceq(2) = x(4) - ((phi.*(v+x(2)+delta.*x(1)) + (1+(x(3)./(1-x(3)))).*(x(1)+x(2)))./Yprod);
c = [];
end

[TerribleStudent]

Sorry for late reply, the plot looks like this for me

[Torsten]

Did you check the exitflags from "fmincon" ? Maybe the peaks correspond to artificial "solutions" where the solver didn't converge.

[TerribleStudent]

No, I did not. I'm not familiar with this procedure/troubleshooting.

[Torsten]

Replace

[Xop, Fop] = fmincon(Objfcn,x0,[],[],[],[],lb,ub,nonlcon,options);

by

[Xop, Fop, exitflag(i,j),~] = fmincon(Objfcn,x0,[],[],[],[],lb,ub,nonlcon,options);

and check the exitflags for the different runs. They are explained in the documentation for "fmincon" and indicate whether a run was successful or not.

[TerribleStudent]

Thank you.