Circular dependency in function and optimization problem

I'm trying to maximize the following optimization problem for "profit" and both inputs

%% Solve Optimization problem

x0 = [2 2];
nonlcon = @constraint;
Objfcn = @optmprob;

options = optimoptions(@fmincon,'Algorithm','interior-point','Display','iter');
[Xop, Fop] = fmincon(Objfcn,x0,[],[],[],[],[],[],nonlcon,options)

where the problem is the circular dependency of some variables in the objective function (all global parameters are defined in the master script):

function profit = optmprob(x)

% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev evP;


Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);

p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+r).*(x(1)+x(2)))./Yprod);

pcdf = cdf("Normal",p_bar,mu,sdev);

r = pcdf./(1-pcdf);


profit = -(pcdf.*phi.*(v+x(2)+delta.*x(1)) + evP.*Yprod - (1-pcdf).*(1+r).*(x(1)+x(2)));

end

p_bar is a variable that depends on both inputs (i.e. x(1) and x(2)) and the rate r, the rate on the otherhand is defiend as a ratio over the cdf of the p_bar variable. Hence, the circular dependency.

I'm getting the error "Unrecognized function or variable 'r'.", and I have tried to circumvent it with handles @ but then I cannot solve it with fmincon due to non double values.

I'm quite new to optimization problems and Matlab in general so all suggestions is highly appreciated regarding how to solve issues like this one and how to proceed forward. Thanks.

Sincerely,

TS

Answer 1

[Torsten]

Define r as a third optimization variable and use the nonlinear constraint function "nonlcon" as

function profit = optmprob(x)
  % Calibration
  global alphaK alphaL bbeta delta phi rho v w mu sdev evP;
  Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
  p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+x(3)).*(x(1)+x(2)))./Yprod);
  pcdf = cdf("Normal",p_bar,mu,sdev);
  r = pcdf./(1-pcdf);
  profit = -(pcdf.*phi.*(v+x(2)+delta.*x(1)) + evP.*Yprod - (1-pcdf).*(1+x(3)).*(x(1)+x(2)));
end
function [c,ceq] = nonlcon(x)
  % Calibration
  global alphaK alphaL bbeta delta phi rho v w mu sdev evP
  Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
  p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+x(3)).*(x(1)+x(2)))./Yprod);
  pcdf = cdf("Normal",p_bar,mu,sdev);
  ceq = x(3) - pcdf./(1-pcdf);
  c = [];
end

[TerribleStudent]

Thank you, that makes perfect sense.

Just to be clear, this implies that I can just replace r with x(3) in my code, which then will look like this giving me a second contraint in addition to my main.

function profit = optmprob(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev evP;

Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+x(3)).*(x(1)+x(2)))./Yprod);
pcdf = cdf("Normal",p_bar,mu,sdev);
%r = pcdf./(1-pcdf);

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


function [c,ceq] = constraint(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev evP;

Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+x(3)).*(x(1)+x(2)))./Yprod);
pcdf = cdf("Normal",p_bar,mu,sdev);
%r = pcdf./(1-pcdf);

ceq(1) = p_bar - ((phi.*(v+x(2)+delta.*x(1)) + (1+x(3)).*(x(1)+x(2)))./Yprod);

ceq(2) = x(3) - pcdf./(1-pcdf);

c = [];

end

However, the starting value for x(3) in x0 is not straightforward, and the following error message imples complex solution for the rate r?

Error using erfc

Input must be real and full.

Error in normcdf>localnormcdf (line 134)

p(todo) = 0.5 * erfc(-z ./ sqrt(2));

Thank you for input.

Best, TS

[Torsten]

It's necessary to have executable code that reproduces the error.

It seems p_bar gets NaN in the course of the calculation so that cdf("Normal",p_bar,mu,sdev); cannot be evaluated.

Further I don't understand why you need to set the first equality constraint. You define

p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+x(3)).*(x(1)+x(2)))./Yprod);

therefore you don't need to force the difference of left and right hand side to be zero. This will always be true.

The below code works technically. Adapt the global variables according to your needs:

%% Solve Optimization problem
global alphaK alphaL bbeta delta phi rho v w mu sdev evP
alphaK =1;
alphaL=1;
bbeta =1;
delta=1;
phi =1;
rho =1;
v=1;
w =1;
mu=1;
sdev=1;
evP=1;
x0 = [0 2 0.5];
nonlcon = @constraint;
Objfcn = @optmprob;
options = optimoptions(@fmincon,'Algorithm','interior-point','Display','iter');
[Xop, Fop] = fmincon(Objfcn,x0,[],[],[],[],[],[],nonlcon,options)

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       4   -5.048068e+00    4.803e+00    1.406e+00
    1       8   -8.096482e+00    1.980e+00    1.513e+00    2.306e+00
    2      12   -4.148399e+02    1.282e+01    2.665e+01    2.296e+01
    3      16   -9.405825e+01    5.135e+00    2.767e+01    2.060e+01
    4      20   -2.268094e+02    2.365e+00    9.955e+00    6.677e+01
    5      24   -7.987852e+02    2.217e+00    3.330e+02    1.715e+02
    6      30   -4.470473e+02    1.924e+00    8.326e+01    7.307e+01
    7      34   -4.085741e+04    1.714e+01    2.378e+03    1.481e+03
    8      38   -2.668906e+02    5.243e+00    3.913e+00    1.574e+03
    9      42   -6.413351e+03    2.169e+00    1.824e+02    2.921e+03
   10      46   -3.556462e+04    2.690e+00    1.386e+04    6.623e+03
   11      51   -1.122996e+04    1.947e+00    1.445e+03    4.422e+03
   12      55   -2.280890e+05    7.009e+00    3.254e+04    1.778e+04
   13      60   -6.196702e+04    3.511e+00    8.846e+02    1.048e+04
   14      65   -1.816043e+04    2.027e+00    2.211e+03    5.341e+03
   15      72   -1.210165e+04    1.930e+00    5.878e+03    1.608e+03
   16     101   -1.249369e+04    1.919e+00    1.581e+00    1.963e-01
Nonlinear constraint function returned Inf; trying a new point...
   17     106   -1.249212e+04    1.919e+00    2.431e+00    3.927e-01
Nonlinear constraint function returned Inf; trying a new point...
   18     118   -1.237888e+04    1.918e+00    1.562e+00    4.908e-02
   19     123   -1.237917e+04    1.918e+00    4.201e+01    9.817e-02
   20     134   -1.243458e+04    1.919e+00    1.574e+00    2.454e-02
Nonlinear constraint function returned Inf; trying a new point...
   21     139   -1.243451e+04    1.919e+00    1.574e+00    2.454e-02
Nonlinear constraint function returned Inf; trying a new point...
   22     148   -1.237916e+04    1.918e+00    1.562e+00    2.454e-02
   23     153   -1.237929e+04    1.918e+00    1.109e+01    4.908e-02
   24     165   -1.239291e+04    1.918e+00    1.571e+00    6.135e-03
Nonlinear constraint function returned Inf; trying a new point...
   25     170   -1.239293e+04    1.918e+00    1.571e+00    1.227e-02
Nonlinear constraint function returned Inf; trying a new point...
   26     181   -1.238611e+04    1.918e+00    1.581e+00    3.068e-03
Nonlinear constraint function returned Inf; trying a new point...
   27     186   -1.238613e+04    1.918e+00    1.582e+00    6.135e-03
Nonlinear constraint function returned Inf; trying a new point...
   28     198   -1.238443e+04    1.918e+00    7.377e+00    7.669e-04
   29     203   -1.238443e+04    1.918e+00    7.377e+00    1.534e-03
   30     222   -1.238443e+04    1.918e+00    3.121e+01    1.498e-06
Nonlinear constraint function returned Inf; trying a new point...

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   31     227   -1.238443e+04    1.918e+00    7.161e+01    2.996e-06
Nonlinear constraint function returned Inf; trying a new point...
   32     239   -1.238443e+04    1.918e+00    2.768e+03    3.745e-07
Nonlinear constraint function returned Inf; trying a new point...
   33     244   -1.238443e+04    1.918e+00    2.760e+03    3.745e-07

Converged to an infeasible point.

fmincon stopped because it is unable to find a point locally that satisfies
the constraints within the value of the constraint tolerance.


Consider enabling the interior point method feasibility mode.

Xop = 1×3

1.0e+03 * 3.9497 3.9517 -0.0013

Fop = -1.2384e+04

function profit = optmprob(x)
  % Calibration
  global alphaK alphaL bbeta delta phi rho v w mu sdev evP;
  Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
  p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+x(3)).*(x(1)+x(2)))./Yprod);
  pcdf = cdf("Normal",p_bar,mu,sdev);
  r = pcdf./(1-pcdf);
  profit = -(pcdf.*phi.*(v+x(2)+delta.*x(1)) + evP.*Yprod - (1-pcdf).*(1+x(3)).*(x(1)+x(2)));
end
function [c,ceq] = constraint(x)
  % Calibration
  global alphaK alphaL bbeta delta phi rho v w mu sdev evP
  Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
  p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+x(3)).*(x(1)+x(2)))./Yprod);
  pcdf = cdf("Normal",p_bar,mu,sdev);
  ceq = x(3) - pcdf./(1-pcdf);
  c = [];
end

[TerribleStudent]

Thank you. I defined p_bar from the constraint, but I was unsure if I still had to specify it in the constraint function for the optimization process to work correctly in Matlab.

My last question is, if I want to multiply one of the terms in the objective function with a random variable, and plot it, how can I do that most effectively?

I have at the moment specified evP as

mu = 1;                                % Mean
sdev = 1;                              % Standard deviation
rng(0,'twister')
evP = sdev.*randn()+mu;

but this only gives me a single scalar, what if I want noise/uncertainty over all possible values i.e. I want several local min/maxima in the 3d plot. Thank you for your time.

[Torsten]

You can't solve problems where inputs change randomly during one run using "fmincon".

"fmincon" is a deterministic solver and not suited for stochastic problems.

Or what do you intend if you say that you "want noise/uncertainty over all possible values" ?

[TerribleStudent]

Yes, so fmincon might not be the suitable solver (maybe patternsearch is better) in this case. I don't want the inputs to change randomly, but rather to make the objetive function subject to some uncertainty (which might imply other optimal points?).

I have now solved the problem determenistic, but I would like to add some uncertainty to the Yprod term in the objective function, similar to this example where they perturb the objectve function Optimization of Stochastic Objective Function - MATLAB & Simulink - MathWorks India

But I just want to add this uncertainty to one of the terms if this is possible, and 3D plot it in a similar way. But I cannot figure it out.

[Torsten]

It's not possible to change Yprod randomly in the course of a computation with "fmincon", as I already explained.

But you could call your problem several times with the parameters that make up Yprod (alphaK,rho , alphaL,v,w,bbeta) changed appropriately.

[TerribleStudent]

Thank you. I have one last question however that I hope you can help to shed some light on. x(4) is defined as the cdf of p_bar in the constraint as well as x(3) = x(4)/(1-x(4)), x(4) should thus be a probability but the optimal value obtained is a really large number, hence no probabiliy. This is my code:

%% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP p_bar;
alphaK = 0.5;
alphaL = 0.5;
bbeta = 1;
delta = 0.5;
phi = 0.8;
rho = 0.2;
v = 0.001;
w = 1;
mu = 1;                                % Mean
sdev = 1;                              % Standard deviation
rng(0,'twister')
p = mu + sdev.*randn(10000,1);
cevP = 1;
p_bar = 1;
%% Solve Optimization problem
x0 = [2 2 2 2];                      % Start values for optimization
nonlcon = @constraint;
Objfcn = @optmprob;
options = optimoptions(@fmincon,'Algorithm','interior-point','Display','iter');
[Xop, Fop] = fmincon(Objfcn,x0,[],[],[],[],[],[],nonlcon,options)


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

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


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

ceq = x(3) - x(4)./(1-x(4));
c = [];
end

When I simply run normcdf(p_bar,mu,sdev) in Matlab window I obtain 0.6378, which is not the optimal value obtained in Xop(4), why?

Thank you very much, it's highly appreciated!

Sincerely.

[Torsten]

x is an input to "constraint". So x(4) has a value from "fmincon" that you can use to set your constraints c and ceq. It's not possible to simply reset its value to normcdf(p_bar,mu,sdev). If you want this relation to hold, you can set a second equality constraints as ceq(2) = x(4) - normcdf(p_bar,mu,sdev).

By the way: The relation ceq = x(3) - x(4)./(1-x(4)) gives x(3) = x(4)./(1-x(4)). Thus in your objective function you could substitute x(3) by x(4)/(1-x(4)) and reduce the number of unknown parameters from 4 to 3.

[TerribleStudent]

When I substitute x(3) for x(4)/(1-x(4)) and reduce the number of variables to 3 (hence replace x(4) with x(3)) I still get different values for Xop and normcdf(p_bar,mu,sdev).

Also there seems to be an issue when I devide the integral in cevP by (1-x(3)) and not (1-normcdf(.)), but it should be the same if I specify x(3) =normcdf() in the constraint? I really don't understand why I get different results.

%% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP p_bar;
alphaK = 0.5;
alphaL = 0.5;
bbeta = 1;
delta = 0.5;
phi = 0.8;
rho = 1;
v = 0.001;
w = 1;
mu = 1;                                % Mean
sdev = 1;                              % Standard deviation
rng(0,'twister')
p = mu + sdev.*randn(10000,1);
cevP = 1;
p_bar = 1;
%% Solve Optimization problem
x0 = [1 1 0.5];                      % Start values for optimization
nonlcon = @constraint;
Objfcn = @optmprob;
options = optimoptions(@fmincon,'Algorithm','interior-point','Display','iter');
[Xop, Fop] = fmincon(Objfcn,x0,[],[],[],[],[],[],nonlcon,options)


function profit = optmprob(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP p_bar;

Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);

p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+(x(3)./(1-x(3)))).*(x(1)+x(2)))./Yprod);

%x(4) = normcdf(p_bar,mu,sdev);

cevP = quadgk(@(p) p.*normpdf(p,mu,sdev)./(1-x(3)), p_bar, Inf);

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


function [c,ceq] = constraint(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev;

Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);

p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+(x(3)./(1-x(3)))).*(x(1)+x(2)))./Yprod);

%x(4) = normcdf(p_bar,mu,sdev);
%ceq(1) = x(3) - (x(3)./(1-x(3)));

ceq = x(3) - normcdf(p_bar,mu,sdev);
c = [];
end

[Torsten]

This would be the code if you replace x(3) by x(4)/(1-x(4)) or - what I did below - x(4) by x(3)/(1+x(3))

Maybe you have to set lower and upper bounds on the variables in order to avoid the upcoming error - I don't know what you are doing.

%% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP p_bar;
alphaK = 0.5;
alphaL = 0.5;
bbeta = 1;
delta = 0.5;
phi = 0.8;
rho = 0.2;
v = 0.001;
w = 1;
mu = 1;                                % Mean
sdev = 1;                              % Standard deviation
rng(0,'twister')
p = mu + sdev.*randn(10000,1);
cevP = 1;
p_bar = 1;
%% Solve Optimization problem
x0 = [2 2 2];                      % Start values for optimization
nonlcon = @constraint;
Objfcn = @optmprob;
options = optimoptions(@fmincon,'Algorithm','interior-point','Display','iter');
[Xop, Fop] = fmincon(Objfcn,x0,[],[],[],[],[],[],nonlcon,options)

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       4   -1.230912e+01    3.333e-01    1.209e+01

Error using erfc
Input must be real and full.

Error in normcdf>iErf (line 150)
p = 0.5 * erfc(-z ./ sqrt(2));

Error in normcdf>localnormcdf (line 128)
p = iErf(z);

Error in normcdf (line 50)
[varargout{1:max(1,nargout)}] = localnormcdf(uflag,x,varargin{:});

Error in solution>constraint (line 37)
ceq = x(3)/(1+x(3)) - normcdf(p_bar,mu,sdev);

Error in barrier

Error in fmincon (line 891)
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

function profit = optmprob(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP p_bar;
Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+x(3)).*(x(1)+x(2)))./Yprod);
%x(4) = normcdf(p_bar,mu,sdev);
cevP = quadgk(@(p) p.*normpdf(p,mu,sdev)./(1-normcdf(p_bar,mu,sdev)), p_bar, Inf);
profit = -(x(3)/(1+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;
Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+x(3)).*(x(1)+x(2)))./Yprod);
ceq = x(3)/(1+x(3)) -  normcdf(p_bar,mu,sdev);
c = [];
end

[TerribleStudent]

Hi again. Another issue that came to light regarding the circular dependency is that I get a NaN value for the variable cevP, which is the conditonal expected value, see below. The problem, as far as I understand it, is related to the relative high value obtained for p_bar, which results in a CDF value of 1. Hence, the denominator in the cevP variable becomes 1-1=0, hence NaN. Here is my code:

%% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP p_bar;
alphaK = 0.5;
alphaL = 0.5;
bbeta = 0.99;
delta = 0.5;
phi = 0.8;
rho = 0.99;
v = 0.001;
w = 1;
mu = 1;                                % Mean
sdev = 1;                              % Standard deviation
rng(0,'twister')
p = mu + sdev.*randn(10000,1);
cevP = 1;
p_bar = 1;
%% Solve Optimization problem
x0 = [1 1 0.5];                          % Start values for optimization
lb = [0 0 0];                          % Lower bound
nonlcon = @constraint;
Objfcn = @optmprob;
options = optimoptions(@fmincon,'Algorithm','interior-point','Display','iter');
[Xop, Fop] = fmincon(Objfcn,x0,[],[],[],[],lb,[],nonlcon,options)

function profit = optmprob(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP p_bar;
Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+(x(3)./(1-x(3)))).*(x(1)+x(2)))./Yprod);

cevP = quadgk(@(p) (p.*normpdf(p,mu,sdev))./(1-normcdf(p_bar,mu,sdev)), p_bar, Inf);

profit = -(x(3).*phi.*(v+x(2)+delta.*x(1)) + cevP.*Yprod - (1-x(3)).*(1+(x(3)./(1-x(3)))).*(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);
p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+(x(3)./(1-x(3)))).*(x(1)+x(2)))./Yprod);

ceq = x(3) - normcdf(p_bar,mu,sdev);
c = [];
end

I think the problem is also related to that x(3) shows up in the numerator of p_bar, hence the optimization tries to maximize it and therefore we receive a high value. But since I have a circular dependency between the CDF ( x(3) ), and p_bar I don't know how to circumvent it. Any ideas? Thanks a lot for your time and effort sir, it's highly appreciated. Have a nice day!

Sincerely, TS

[Torsten]

Why do your equations change every time you post a new code suggestion ?

I don't know the underlying problem. So I think it doesn't make sense to search for something from which I don't know what it is.

If you can turn cevP as large as you want by making 1-normcdf(p_bar,mu,sdev) zero or almost zero, maybe your problem is unbounded.

By the way: The expression (1-x(3)).*(1+(x(3)./(1-x(3)))) in profit equals 1.

[TerribleStudent]

I edit the equations since I have made some simplificatons or redefined the definition of some of the variables (and I still havent solved it). My issue is that I do not know how to specify the variables in an optimization problem when some variables, that are not part of the optimization problem, is a function of each other, i.e. r is a function of F(p), p contains r, and F(p) is the cdf of p. I don't want the algorithm to optimize these variables per se, they should rather be seen as exogenous and still hold for some constraint.

Now I have a variable cevP defined as conditional expectation, in integral form, but it still doesn't equal, or larger to, its condition.

[Torsten]

My issue is that I do not know how to specify the variables in an optimization problem when some variables, that are not part of the optimization problem, is a function of each other, i.e. r is a function of F(p), p contains r, and F(p) is the cdf of p. I don't want the algorithm to optimize these variables per se, they should rather be seen as exogenous and still hold for some constraint.

All variables that are part of your optimization, that can't be written as a function of the optimization variables and for which you want to prescribe constraints have to be defined as optimization variables as well. If the solver doesn't respect your constraint, you should try to start from an initial value for x0(3) that satisfies your constraint. Then you should check why your objective function returns NaN for the x0 vector that satisfies your constraint (see below).

%% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP p_bar;
alphaK = 0.5;
alphaL = 0.5;
bbeta = 0.99;
delta = 0.5;
phi = 0.8;
rho = 0.99;
v = 0.001;
w = 1;
mu = 1;                                % Mean
sdev = 1;                              % Standard deviation
rng(0,'twister')
p = mu + sdev.*randn(10000,1);
cevP = 1;
p_bar = 1;
%% Solve Optimization problem
x0 = [1 2 0.5]; 
for i = 1:20
Yprod = (alphaK.*x0(1).^rho + alphaL.*((v + x0(2)) ./w).^rho).^(bbeta./rho);
p_bar = ((phi.*(v+x0(2)+delta.*x0(1)) + (x0(1)+x0(2)))/(1-x0(3))./Yprod);
x0(3) = normcdf(p_bar,mu,sdev);
end
x0(3)

ans = 1

Yprod = (alphaK.*x0(1).^rho + alphaL.*((v + x0(2)) ./w).^rho).^(bbeta./rho);
p_bar = ((phi.*(v+x0(2)+delta.*x0(1)) + (x0(1)+x0(2)))/(1-x0(3))./Yprod);
normcdf(p_bar,mu,sdev)

ans = 1

optmprob(x0)

Warning: Infinite or Not-a-Number value encountered.

ans = NaN

lb = [0 0 0];                          % Lower bound
nonlcon = @constraint;
Objfcn = @optmprob;
options = optimoptions(@fmincon,'Algorithm','interior-point','Display','iter');
[Xop, Fop] = fmincon(Objfcn,x0,[],[],[],[],lb,[],nonlcon,options)

Warning: Infinite or Not-a-Number value encountered.

Error using barrier
Objective function is undefined at initial point. Fmincon cannot continue.

Error in fmincon (line 891)
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

function profit = optmprob(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP p_bar;
Yprod = (alphaK.*x(1).^rho + alphaL.*((v + x(2)) ./w).^rho).^(bbeta./rho);
p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+(x(3)./(1-x(3)))).*(x(1)+x(2)))./Yprod);
cevP = quadgk(@(p) (p.*normpdf(p,mu,sdev))./(1-normcdf(p_bar,mu,sdev)), p_bar, Inf);
profit = -(x(3).*phi.*(v+x(2)+delta.*x(1)) + cevP.*Yprod - (1-x(3)).*(1+(x(3)./(1-x(3)))).*(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);
p_bar = ((phi.*(v+x(2)+delta.*x(1)) + (1+(x(3)./(1-x(3)))).*(x(1)+x(2)))./Yprod);
ceq = x(3) - normcdf(p_bar,mu,sdev);
c = [];
end

[TerribleStudent]

All variables that are part of your optimization, that can't be written as a function of the optimization variables and for which you want to prescribe constraints have to be defined as optimization variables as well.

Thank you, this makes sense. I think I might have solved it. So in my case, I have two variables that I want to optimize, x1 and x2, but since I also have the CDF and p_bar that has to fullfill some constraint, I can declare these as optimizations variables x3 and x4, respectively, as well. The constraint will thus look like

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);

and then the cevP variable can be defined as

cevP = quadgk(@(p) (p.*normpdf(p,mu,sdev))./(1-x(3)), x(4), Inf);

In this case, I don't have to declare cevP as an additional variable,x5,? This gives reasonable result except that the resulting optimal values for x1 and x2 are, or close to, 0. Thanks.

[Torsten]

Usually, it's not necessary to declare x(4) or cevP as optimization variables because they can be explicitly expressed by x(1)-x(3). But if your code works better this way ...

[TerribleStudent]

The code works when x4 is declared as above, in contrast to just x3.

My last though is that, however, when I plot the function and optimal solution, this seems not to be the maximum point even though the problem is defined as to minimize the negative function:

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

why is that?

[Torsten]

My last though is that, however, when I plot the function and optimal solution, this seems not to be the maximum point

Why do you think it's not the optimum ?

What is your actual code ?

How do you plot the function ? It depends on three inputs x(1)-x(3) !

[TerribleStudent]

Why do you think it's not the optimum ?

Because the surface plots indicate that it's a (local) minimum. However, given the calibration/ constraints it might be a maximum I guess?

What is your actual code ?

How do you plot the function ? It depends on three inputs x(1)-x(3) !

I plot the objective function on the z-axis given the obtained optimal value for x3, and x1 and x2 are x and y-axis, respectively. Then I mark the optimal point given x1, x2 and x3.

[X,Y] = meshgrid(0: 0.1: 10);          % Range to plot
figure;
Z = Xop(3).*phi.*(v+Y+delta.*X) + cevP.*(alphaK.*X.^rho + alphaL.*((v + Y) ./w).^rho).^(bbeta./rho) - (1-Xop(3)).*(1+(Xop(3)./(1-Xop(3)))).*(X+Y);
surf(X,Y,Z)
xlabel('K');
ylabel('B');
zlabel('');
hold on
scatter3(Xop(1),Xop(2),-Fop,150,'filled','o','r')
text(Xop(1),Xop(2),Fop,'Optimum', 'FontSize',12)
hold off

According to the plot, I higher value can be optained by increasing x1 and x2, however, I guess the constraints is restrcting it??

[Torsten]

You mean a plot like this ?

%% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP p_bar;
alphaK = 0.5;
alphaL = 0.5;
bbeta = 0.99;
delta = 0.5;
phi = 0.8;
rho = 0.99;
v = 0.001;
w = 1;
mu = 1;                                % Mean
sdev = 1;                              % Standard deviation
rng(0,'twister')
p = mu + sdev.*randn(10000,1);
cevP = 1;
p_bar = 1;
%% Solve Optimization problem
x0 = [1 1 0.5 0.5]; 
lb = [0 0 0 0];                          % Lower bound
nonlcon = @constraint;
Objfcn = @optmprob;
options = optimoptions(@fmincon,'Algorithm','interior-point','Display','iter');
[Xop, Fop] = fmincon(Objfcn,x0,[],[],[],[],lb,[],nonlcon,options)

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       5   -6.884887e-01    4.698e+00    3.784e-01
    1      11   -2.072775e-01    4.176e+00    1.128e+00    2.930e-01
    2      16    1.851202e-01    3.745e+00    7.698e-01    8.834e-01
    3      22    2.086802e-01    3.682e+00    7.164e-01    1.914e-01
    4      27    2.271461e-01    3.622e+00    5.899e-01    1.625e-01
    5      33    8.246304e-02    3.629e+00    3.449e-01    2.426e-01
    6      38    1.973537e-01    3.554e+00    4.928e-01    3.697e-01
    7      43   -2.875754e-04    3.130e+00    9.868e+01    1.433e+00
    8      48   -1.012412e-03    1.512e+00    2.032e+03    2.049e-02
    9      53   -1.010477e-03    1.442e+00    2.629e+02    9.580e-03
   10      58   -1.010087e-03    1.436e+00    2.819e+02    8.883e-04
   11      64   -1.010142e-03    1.436e+00    2.827e+02    2.018e-04
   12      69   -1.206631e-03    1.487e-01    3.146e+03    1.461e+00
   13      74   -1.696113e-03    7.411e-05    6.843e+03    1.488e-01
   14      79   -1.696584e-03    3.796e-07    6.848e+03    6.475e-04
   15      84   -1.697597e-03    2.130e-06    6.860e+03    1.373e-03
   16      89   -1.699792e-03    9.672e-06    6.886e+03    2.979e-03
   17      94   -1.704900e-03    4.889e-05    6.946e+03    6.929e-03
   18      99   -1.716288e-03    2.349e-04    7.083e+03    1.541e-02
   19     105   -1.741806e-03    4.295e-05    7.417e+03    3.517e-02
   20     111   -1.782764e-03    2.823e-05    8.018e+03    5.697e-02
   21     117   -1.793146e-03    2.772e-05    8.186e+03    1.440e-02
   22     123   -1.821591e-03    7.805e-05    8.674e+03    3.917e-02
   23     129   -1.888517e-03    8.752e-05    1.004e+04    9.418e-02
   24     135   -1.901465e-03    2.527e-05    1.035e+04    1.844e-02
   25     144   -1.906275e-03    3.083e-05    1.047e+04    6.703e-03
   26     157   -1.917198e-03    7.443e-06    8.622e+00    1.569e-02
   27     174   -1.927482e-03    6.743e-06    8.767e+00    1.481e-02
   28     186   -1.938519e-03    8.692e-06    8.928e+00    1.596e-02
   29     198   -1.950403e-03    1.120e-05    7.734e-02    1.727e-02
   30     204   -2.480830e-03    3.953e-01    9.395e+01    3.845e-01

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   31     209   -2.157094e-03    1.351e-01    4.424e+01    1.079e-01
   32     214   -2.082969e-03    3.233e-02    2.622e+01    9.905e-02
   33     219   -2.054003e-03    2.452e-03    2.157e+01    2.570e-02
   34     224   -2.058832e-03    6.457e-05    2.301e+01    8.449e-03
   35     229   -2.086717e-03    2.442e-03    3.030e+01    4.013e-02
   36     235   -2.190281e-03    8.092e-04    6.642e+01    1.562e-01
   37     240   -2.195318e-03    6.625e-04    2.199e+01    2.134e-02
   38     245   -2.198060e-03    4.508e-05    8.027e+00    5.751e-03
   39     250   -2.198188e-03    1.136e-07    1.698e-01    3.019e-04
   40     255   -2.198181e-03    2.231e-10    2.114e-02    1.060e-05
   41     260   -2.200548e-03    2.525e-05    4.260e+00    3.939e-03
   42     265   -2.210496e-03    4.312e-04    6.636e-01    1.611e-02
   43     270   -2.209812e-03    7.998e-07    4.962e-01    4.197e-04
   44     275   -2.209908e-03    4.011e-08    6.532e-02    1.608e-04
   45     280   -2.210003e-03    4.086e-08    5.284e-03    1.598e-04
   46     285   -2.210007e-03    6.189e-11    4.000e-03    7.420e-06
   47     290   -2.221857e-03    6.164e-04    5.502e+00    1.910e-02
   48     295   -2.282556e-03    1.299e-02    4.669e-01    8.277e-02
   49     300   -2.261089e-03    8.876e-04    1.378e+00    1.451e-02
   50     305   -2.264675e-03    6.886e-05    8.184e-01    6.884e-03
   51     310   -2.270373e-03    1.678e-04    2.323e-01    9.882e-03
   52     315   -2.270533e-03    4.147e-07    2.840e-02    6.173e-04
   53     320   -2.270529e-03    7.061e-11    8.000e-04    4.424e-06
   54     325   -2.322098e-03    1.004e-02    5.704e+00    7.177e-02
   55     330   -2.754478e-03    2.429e-01    4.545e+00    2.707e-01
   56     335   -2.437336e-03    6.560e-02    5.982e+00    5.016e-02
   57     340   -2.438020e-03    1.827e-03    4.120e+00    1.858e-02
   58     345   -2.440522e-03    8.697e-05    2.519e+00    6.108e-03
   59     350   -2.453291e-03    9.965e-04    2.093e+00    2.225e-02
   60     355   -2.493168e-03    7.669e-03    7.982e-01    5.906e-02

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   61     360   -2.497484e-03    1.013e-03    1.977e-01    2.736e-02
   62     367   -2.510653e-03    1.782e-03    7.360e-01    2.472e-02
   63     374   -2.526931e-03    2.758e-03    9.262e-01    2.978e-02
   64     379   -2.547879e-03    3.348e-03    3.171e-01    4.036e-02
   65     384   -2.555667e-03    1.043e-03    7.121e-01    2.414e-02
   66     391   -2.563690e-03    1.098e-03    1.540e-01    1.652e-02
   67     396   -2.578232e-03    1.695e-03    8.710e-01    2.800e-02
   68     401   -2.579659e-03    1.293e-04    9.778e-01    8.922e-03
   69     407   -2.585488e-03    4.653e-05    1.831e-01    1.285e-02
   70     412   -2.591465e-03    2.836e-04    5.704e-03    1.119e-02
   71     417   -2.591500e-03    1.863e-06    4.740e-02    1.144e-03
   72     423   -2.593513e-03    3.355e-06    9.987e-02    4.200e-03
   73     429   -2.595304e-03    3.342e-06    5.197e-02    3.718e-03
   74     434   -2.596461e-03    1.215e-05    1.440e-01    2.342e-03
   75     439   -2.596830e-03    1.385e-06    1.112e-01    8.079e-04
   76     445   -2.597328e-03    9.815e-08    6.861e-02    1.040e-03
   77     450   -2.597604e-03    7.043e-07    4.514e-02    5.715e-04
   78     455   -2.597712e-03    1.027e-07    3.593e-02    2.246e-04
   79     460   -2.597796e-03    6.163e-08    2.873e-02    1.752e-04
   80     465   -2.597837e-03    1.287e-08    2.523e-02    8.523e-05
   81     470   -2.597856e-03    2.234e-09    2.353e-02    4.139e-05
   82     475   -2.597866e-03    1.124e-10    2.273e-02    1.964e-05
   83     481   -2.597868e-03    7.024e-11    2.254e-02    4.548e-06
   84     494   -2.597868e-03    1.552e-11    7.741e-05    4.187e-07
   85     499   -1.238825e-02    5.699e+00    1.947e+02    4.012e-01
   86     504   -5.195370e-03    5.571e+00    7.540e+00    1.017e+00
   87     509   -4.982105e-03    5.199e+00    7.356e+00    2.120e-02
   88     514   -2.808766e-03    4.950e-01    5.591e+00    2.887e-01
   89     519   -2.797781e-03    4.682e-01    5.688e+00    3.996e-03
   90     524   -2.580018e-03    5.538e-02    4.000e+00    1.301e-01

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   91     529   -2.586742e-03    2.532e-02    2.299e+00    8.469e-02
   92     534   -2.583970e-03    1.368e-03    1.300e+00    4.216e-02
   93     539   -2.590120e-03    1.210e-03    6.452e-01    2.383e-02
   94     544   -2.595823e-03    5.172e-04    1.870e-01    1.541e-02
   95     552   -2.599194e-03    5.131e-04    4.277e-02    7.326e-03
   96     564   -2.611032e-03    9.243e-05    5.368e-05    2.653e-02
   97     569   -2.611081e-03    3.267e-07    1.264e-03    5.012e-05
   98     574   -2.612606e-03    8.516e-06    8.238e-02    2.292e-03
   99     585   -2.620144e-03    2.581e-05    5.199e-05    1.597e-02
  100     594   -2.622627e-03    3.343e-05    1.785e-01    4.955e-03
  101     605   -2.627279e-03    5.796e-06    5.080e-05    9.967e-03
  102     619   -2.635592e-03    3.234e-05    4.948e-05    1.735e-02
  103     631   -2.643802e-03    3.156e-05    4.822e-05    1.719e-02
  104     644   -2.651923e-03    3.068e-05    4.703e-05    1.702e-02
  105     657   -2.659959e-03    2.985e-05    4.590e-05    1.685e-02
  106     670   -2.667912e-03    2.907e-05    4.482e-05    1.669e-02
  107     683   -2.675787e-03    2.834e-05    4.380e-05    1.654e-02
  108     695   -2.683587e-03    2.764e-05    4.282e-05    1.639e-02
  109     707   -2.691314e-03    2.699e-05    4.189e-05    1.624e-02
  110     720   -2.698972e-03    2.637e-05    4.100e-05    1.610e-02
  111     733   -2.706562e-03    2.578e-05    4.015e-05    1.596e-02
  112     746   -2.714088e-03    2.522e-05    3.934e-05    1.583e-02
  113     759   -2.721552e-03    2.468e-05    3.856e-05    1.570e-02
  114     772   -2.728955e-03    2.418e-05    3.781e-05    1.557e-02
  115     785   -2.736301e-03    2.369e-05    3.709e-05    1.545e-02
  116     798   -2.743591e-03    2.323e-05    3.640e-05    1.533e-02
  117     811   -2.750826e-03    2.279e-05    3.573e-05    1.521e-02
  118     824   -2.758009e-03    2.237e-05    3.509e-05    1.510e-02
  119     837   -2.765140e-03    2.197e-05    3.448e-05    1.499e-02
  120     850   -2.772223e-03    2.158e-05    3.388e-05    1.488e-02

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
  121     863   -2.779258e-03    2.121e-05    3.331e-05    1.478e-02
  122     876   -2.786246e-03    2.086e-05    3.276e-05    1.467e-02
  123     889   -2.793189e-03    2.052e-05    3.222e-05    1.457e-02
  124     902   -2.800088e-03    2.019e-05    3.171e-05    1.448e-02
  125     915   -2.806945e-03    1.987e-05    3.130e-05    1.438e-02
  126     928   -2.813760e-03    1.957e-05    3.133e-05    1.429e-02
  127     941   -2.820534e-03    1.928e-05    3.136e-05    1.419e-02
  128     954   -2.827270e-03    1.899e-05    3.139e-05    1.410e-02
  129     967   -2.833966e-03    1.872e-05    3.142e-05    1.402e-02
  130     980   -2.840625e-03    1.846e-05    3.145e-05    1.393e-02
  131     993   -2.847248e-03    1.820e-05    3.148e-05    1.384e-02
  132    1006   -2.853835e-03    1.796e-05    3.150e-05    1.376e-02
  133    1019   -2.860387e-03    1.772e-05    3.153e-05    1.368e-02
  134    1032   -2.866905e-03    1.749e-05    3.155e-05    1.360e-02
  135    1045   -2.873389e-03    1.726e-05    3.157e-05    1.352e-02
  136    1058   -2.879841e-03    1.705e-05    3.159e-05    1.344e-02
  137    1071   -2.886261e-03    1.684e-05    3.161e-05    1.337e-02
  138    1084   -2.892649e-03    1.663e-05    3.163e-05    1.329e-02
  139    1097   -2.899007e-03    1.643e-05    3.164e-05    1.322e-02
  140    1110   -2.905334e-03    1.624e-05    3.166e-05    1.315e-02
  141    1123   -2.911633e-03    1.605e-05    3.168e-05    1.308e-02
  142    1136   -2.917902e-03    1.587e-05    3.169e-05    1.301e-02
  143    1149   -2.924142e-03    1.569e-05    3.170e-05    1.294e-02
  144    1162   -2.930355e-03    1.552e-05    3.172e-05    1.287e-02
  145    1175   -2.936541e-03    1.535e-05    3.173e-05    1.281e-02
  146    1180   -3.160529e-03    6.081e-02    2.308e+01    1.169e-01
  147    1185   -2.991300e-03    2.302e-02    6.569e+00    5.387e-02
  148    1190   -2.961645e-03    2.102e-03    6.083e+00    1.534e-02
  149    1195   -2.959252e-03    8.348e-06    6.059e+00    6.877e-04
  150    1200   -2.961546e-03    5.697e-05    6.176e+00    4.388e-03

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
  151    1205   -2.974386e-03    1.187e-03    6.236e+00    1.950e-02
  152    1212   -2.994216e-03    3.097e-03    3.177e+00    2.818e-02
  153    1219   -3.016889e-03    4.912e-03    5.367e-01    3.286e-02
  154    1224   -3.070697e-03    1.185e-02    1.364e+00    6.008e-02
  155    1229   -3.040807e-03    1.635e-04    1.141e+00    1.663e-02
  156    1234   -3.070488e-03    4.550e-03    4.641e-01    3.626e-02
  157    1239   -3.123733e-03    1.103e-02    8.167e-01    5.687e-02
  158    1244   -3.089999e-03    1.133e-05    9.836e-01    1.056e-02
  159    1249   -3.111014e-03    2.595e-03    1.694e-01    2.720e-02
  160    1254   -3.167029e-03    1.055e-02    9.801e-01    5.365e-02
  161    1259   -3.131361e-03    3.278e-06    9.389e-01    5.027e-03
  162    1264   -3.147610e-03    1.647e-03    2.215e-01    2.155e-02
  163    1269   -3.199791e-03    9.088e-03    1.134e+00    4.899e-02
  164    1274   -3.169440e-03    2.121e-05    1.027e+00    6.057e-03
  165    1279   -3.187542e-03    1.929e-03    1.064e-01    2.298e-02
  166    1284   -3.235043e-03    8.090e-03    1.516e+00    4.615e-02
  167    1289   -3.208787e-03    1.113e-05    1.235e+00    8.852e-03
  168    1294   -3.230935e-03    2.615e-03    2.168e-01    2.631e-02
  169    1299   -3.271929e-03    7.017e-03    2.084e+00    4.327e-02
  170    1304   -3.252930e-03    3.159e-04    1.339e+00    1.467e-02
  171    1309   -3.287489e-03    4.952e-03    1.177e+00    3.534e-02
  172    1314   -3.301531e-03    3.445e-03    3.714e+00    3.284e-02
  173    1319   -3.344859e-03    8.226e-03    9.249e-02    4.630e-02
  174    1324   -3.314645e-03    7.254e-05    8.204e-01    1.112e-02
  175    1329   -3.334959e-03    2.336e-03    1.464e+00    2.431e-02
  176    1334   -3.383236e-03    7.884e-03    2.785e+00    4.407e-02
  177    1339   -3.349475e-03    8.128e-06    2.680e+00    4.288e-03
  178    1344   -3.363542e-03    1.227e-03    3.347e-01    1.763e-02
  179    1349   -3.408190e-03    6.558e-03    3.459e+00    3.956e-02
  180    1354   -3.381235e-03    1.343e-05    3.053e+00    5.421e-03

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
  181    1359   -3.397898e-03    1.568e-03    1.691e-01    1.970e-02
  182    1364   -3.436791e-03    5.664e-03    5.045e+00    3.698e-02
  183    1369   -3.415994e-03    1.020e-04    3.628e+00    9.686e-03
  184    1374   -3.441226e-03    2.882e-03    2.086e+00    2.620e-02
  185    1379   -3.466584e-03    4.149e-03    7.943e+00    3.282e-02
  186    1384   -3.473518e-03    2.502e-03    1.134e+00    2.728e-02
  187    1389   -3.507548e-03    5.671e-03    1.502e+00    3.724e-02
  188    1394   -3.488849e-03    4.255e-04    6.002e+00    1.474e-02
  189    1399   -3.518068e-03    3.655e-03    6.937e+00    2.915e-02
  190    1404   -3.527053e-03    2.300e-03    1.650e+01    2.572e-02
  191    1409   -3.563492e-03    5.733e-03    1.437e+00    3.691e-02
  192    1414   -3.538315e-03    1.040e-04    5.473e+00    1.009e-02
  193    1419   -3.559077e-03    2.119e-03    8.443e+00    2.232e-02
  194    1424   -3.593661e-03    4.901e-03    1.363e+01    3.410e-02
  195    1429   -3.572146e-03    1.214e-04    8.684e+00    9.942e-03
  196    1434   -3.595704e-03    2.410e-03    8.332e+00    2.367e-02
  197    1439   -3.619019e-03    3.442e-03    2.211e+01    2.941e-02
  198    1444   -3.622856e-03    1.846e-03    2.987e+00    2.342e-02
  199    1449   -3.654507e-03    4.534e-03    5.704e+00    3.280e-02
  200    1454   -3.636203e-03    2.798e-04    1.827e+01    1.269e-02
  201    1459   -3.663276e-03    2.859e-03    1.942e+01    2.583e-02
  202    1464   -3.670602e-03    1.725e-03    4.685e+01    2.255e-02
  203    1469   -3.704680e-03    4.558e-03    4.061e+00    3.292e-02
  204    1474   -3.680089e-03    2.218e-05    1.498e+01    8.643e-03
  205    1479   -3.698957e-03    1.514e-03    2.246e+01    1.966e-02
  206    1484   -3.731131e-03    3.784e-03    3.817e+01    3.046e-02
  207    1489   -3.710506e-03    2.477e-05    2.411e+01    8.889e-03
  208    1494   -3.732902e-03    1.783e-03    2.443e+01    2.153e-02
  209    1499   -3.753037e-03    2.416e-03    6.384e+01    2.598e-02
  210    1504   -3.760412e-03    1.595e-03    7.888e+00    2.295e-02

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
  211    1509   -3.787589e-03    3.238e-03    3.464e+00    2.971e-02
  212    1514   -3.771679e-03    1.453e-04    3.920e+01    1.315e-02
  213    1519   -3.798177e-03    2.067e-03    5.312e+01    2.452e-02
  214    1524   -3.800371e-03    7.452e-04    1.380e+02    1.909e-02
  215    1529   -3.835912e-03    3.239e-03    3.358e+01    3.058e-02
  216    1534   -3.809702e-03    1.631e-04    5.986e+01    7.011e-03
  217    1539   -3.826338e-03    5.061e-04    4.755e+01    1.702e-02
  218    1544   -3.858968e-03    2.292e-03    1.071e+02    2.829e-02
  219    1549   -3.836670e-03    2.300e-04    7.515e+01    7.015e-03
  220    1554   -3.856080e-03    3.940e-04    5.005e+01    1.846e-02
  221    1559   -3.880607e-03    1.213e-03    1.656e+02    2.553e-02
  222    1564   -3.873795e-03    2.760e-04    5.299e+01    1.540e-02
  223    1569   -3.906121e-03    1.130e-03    1.059e+02    2.706e-02
  224    1574   -3.890590e-03    5.916e-04    2.197e+02    1.187e-02
  225    1579   -3.923885e-03    4.043e-04    1.679e+02    2.593e-02
  226    1586   -3.920288e-03    6.402e-05    1.534e+01    3.593e-03
  227    1594   -3.958868e-03    1.444e-04    2.129e+02    2.686e-02
  228    1601   -3.947760e-03    3.350e-06    1.226e+02    4.499e-05
  229    1617   -3.946944e-03    1.294e-06    1.655e-05    6.456e-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.

Xop = 1×4

0.0000 0.0000 1.0000 5.6139

Fop = -0.0039

x = 0:0.1:10;
y = x;
[X,Y] = meshgrid(x);          % Range to plot
figure;
max_Z = -Inf;
for i = 1:numel(x)
  for j = 1:numel(y)
    Yprod = (alphaK.*x(i).^rho + alphaL.*((v + y(j)) ./w).^rho).^(bbeta./rho);
    P_bar = ((phi.*(v+y(j)+delta.*x(i)) + (1+(Xop(3)./(1-Xop(3)))).*(x(i)+y(j)))./Yprod);
    cevP = quadgk(@(p) (p.*normpdf(p,mu,sdev))./(1-Xop(3)), P_bar, Inf);
    Z(j,i) = Xop(3).*phi.*(v+y(j)+delta.*x(i)) + cevP.*(alphaK.*x(i).^rho + alphaL.*((v + y(j)) ./w).^rho).^(bbeta./rho) - (1-Xop(3)).*(1+(Xop(3)./(1-Xop(3)))).*(x(i)+y(j));
    if Z(j,i) > max_Z
      max_Z = Z(j,i);
      imax = i;
      jmax = j;
    end
  end
end
x(imax)

ans = 0

y(jmax)

ans = 0

Yprod = (alphaK.*x(imax).^rho + alphaL.*((v + y(jmax)) ./w).^rho).^(bbeta./rho);
P_bar = ((phi.*(v+y(jmax)+delta.*x(imax)) + (1+(Xop(3)./(1-Xop(3)))).*(x(imax)+y(jmax)))./Yprod);
Z(jmax,imax)

ans = 181.9107

format long
[x(imax) y(jmax) Xop(3) P_bar]

ans = 1×4

0 0 0.999998043837827 1.493206881275186

Xop

Xop = 1×4

0.000000002159294 0.000000002159302 0.999998043837827 5.613856436842379

-optmprob([x(imax) y(jmax) Xop(3) P_bar])

ans = 
     1.819106943507891e+02

-optmprob(Xop)

ans = 
   0.003946944141375

surf(X,Y,Z)
xlabel('K');
ylabel('B');
zlabel('');
hold on
scatter3(Xop(1),Xop(2),-Fop,150,'filled','o','r')
text(Xop(1),Xop(2),Fop,'Optimum', 'FontSize',12)
hold off

function profit = optmprob(x)
% Calibration
global alphaK alphaL bbeta delta phi rho v w mu sdev cevP p_bar;
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);
profit = -(x(3).*phi.*(v+x(2)+delta.*x(1)) + cevP.*Yprod - (1-x(3)).*(1+(x(3)./(1-x(3)))).*(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]

It's working now. Thank you for your help, it's highly appreciated.