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
采纳的回答
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)
11 个评论
TerribleStudent
2024-2-21
编辑：TerribleStudent
2024-2-21
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
2024-2-23
Torsten
2024-2-23
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
2024-3-6
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
2024-4-12
Sorry for late reply, the plot looks like this for me 
Torsten
2024-4-12
Did you check the exitflags from "fmincon" ? Maybe the peaks correspond to artificial "solutions" where the solver didn't converge.
TerribleStudent
2024-4-15
Torsten
2024-4-15
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
2024-4-16
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Surrogate Optimization 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
