Cannot get fmincon to work

Hello, I am trying to find the max value of a function of 5 variables with all of them constrained between 2 values. However, I do not have experience using this function and I can't find any good information to use it properly in this scenario either in the MatLab help page, nor in any online tutorial.

As you can see in the code, the function I want to optimize is P, and the variables are V, na, nd, wn and wp. (And yes, I am aware that I may have overused the '.' operator, I was getting the same error over and over again and I decided to nuke it since I don't really care about resources in this case.)

Apologies if the code is too messy to understand.

clear
close
clc
format long

%% Variables
dp = 11.6;
taup = 3710e-6;
dn = 2;
taun = 371e-6;
Jl = 50e-3;
ni = 9.696e9;
q = 1.602176634e-19;
k = 1.3880649e-23;
T = 300;
ep = 1.035918e-12;

vint = @(na,nd) k.*T/q.*log((na.*nd)./ni^2);
xn = @(na,nd) sqrt((2.*ep.*vint(na,nd).*na)./(q.*nd.*(na+nd)));
xp = @(na,nd) sqrt((2.*ep.*vint(na,nd).*nd)./(q.*na.*(na+nd)));

dnxp = @(V,na) ni^2./na.*(exp(q*V./(k*T))-1);
dpxn = @(V,nd) ni^2./nd.*(exp(q*V./(k*T))-1);

%% Functions
C1 = @(V,na,nd,wn) dnxp(V,na)./(-exp((-2.*wn-xp(na,nd))./sqrt(dn*taun))+exp(xp(na,nd)./sqrt(dn*taun)));
C3 = @(V,na,nd,wp) dpxn(V,nd)./(exp(-xn(na,nd)./sqrt(taup*dp))-exp((2.*wp+xn(na,nd))./sqrt(taup*dp)));
Jn = @(x,V,na,nd,wn) q*sqrt(dn/taun).*C1(V,na,nd,wn).*(exp((-2.*wn-x)./sqrt(dn*taun))+exp(x./sqrt(dn*taun)));
Jp = @(x,V,na,nd,wp) -q*sqrt(dp/taup).*C3(V,na,nd,wp).*(exp(x./sqrt(taup*dp))+exp((2.*wp-x)./sqrt(taup*dp)));
Jd = @(V,na,nd,wn,wp) Jp(-xn(na,nd),V,na,nd,wn)+Jn(xp(na,nd),V,na,nd,wp);
J = @(V,na,nd,wn,wp) Jl-Jd(V,na,nd,wn,wp);
P = @(V,na,nd,wn,wp) J(V,na,nd,wn,wp).*V;

%% Boundaries
wn = [0.05e-4, 0.75e-4];
wp = [100e-4, 200e-4];
nd = [1e18, 5e20];
na = [1e14, 1e16];
V = [0, 0.52];

P0 = [(V(1)+V(2))/2; (wn(1)+wn(2))/2; (wp(1)+wp(2))/2; (nd(1)+nd(2))/2; (na(1)+na(2))/2];

%% Optimization
mx = fmincon(-P(V,wn,wp,nd,na),P0,[],[],[],[],[min(V),wn(1),wp(1),nd(1),na(1)],[max(V),wn(2),wp(2),nd(2),na(2)]);

Answer 1

[Torsten]

clear
close
clc
format long
%% Variables
dp = 11.6;
taup = 3710e-6;
dn = 2;
taun = 371e-6;
Jl = 50e-3;
ni = 9.696e9;
q = 1.602176634e-19;
k = 1.3880649e-23;
T = 300;
ep = 1.035918e-12;
vint = @(na,nd) k.*T/q.*log((na.*nd)./ni^2);
xn = @(na,nd) sqrt((2.*ep.*vint(na,nd).*na)./(q.*nd.*(na+nd)));
xp = @(na,nd) sqrt((2.*ep.*vint(na,nd).*nd)./(q.*na.*(na+nd)));
dnxp = @(V,na) ni^2./na.*(exp(q*V./(k*T))-1);
dpxn = @(V,nd) ni^2./nd.*(exp(q*V./(k*T))-1);
%% Functions
C1 = @(V,na,nd,wn) dnxp(V,na)./(-exp((-2.*wn-xp(na,nd))./sqrt(dn*taun))+exp(xp(na,nd)./sqrt(dn*taun)));
C3 = @(V,na,nd,wp) dpxn(V,nd)./(exp(-xn(na,nd)./sqrt(taup*dp))-exp((2.*wp+xn(na,nd))./sqrt(taup*dp)));
Jn = @(x,V,na,nd,wn) q*sqrt(dn/taun).*C1(V,na,nd,wn).*(exp((-2.*wn-x)./sqrt(dn*taun))+exp(x./sqrt(dn*taun)));
Jp = @(x,V,na,nd,wp) -q*sqrt(dp/taup).*C3(V,na,nd,wp).*(exp(x./sqrt(taup*dp))+exp((2.*wp-x)./sqrt(taup*dp)));
Jd = @(V,na,nd,wn,wp) Jp(-xn(na,nd),V,na,nd,wn)+Jn(xp(na,nd),V,na,nd,wp);
J = @(V,na,nd,wn,wp) Jl-Jd(V,na,nd,wn,wp);
P = @(V,na,nd,wn,wp) J(V,na,nd,wn,wp).*V;
%% Boundaries
wn = [0.05e-4, 0.75e-4];
wp = [100e-4, 200e-4];
nd = [1e18, 5e20];
na = [1e14, 1e16];
V = [0, 0.52];
% Either use the order of variables like in the vector of initial
% values and lower and upper bounds in the call to PP or 
% change the order of the variables in the initial value vector P0 and 
% in the lower and upper bounds to (V,na,nd,wn,wp) instead of (V,wn,wp,nd,na).
% I did the first, but I recommend the last (commented out) for clarity reasons.
P0 = [(V(1)+V(2))/2; (wn(1)+wn(2))/2; (wp(1)+wp(2))/2; (nd(1)+nd(2))/2; (na(1)+na(2))/2];
PP = @(p)-P(p(1),p(4),p(5),p(3),p(2));
lb = [min(V),wn(1),wp(1),nd(1),na(1)];
ub = [max(V),wn(2),wp(2),nd(2),na(2)];
%P0 = [(V(1)+V(2))/2; (na(1)+na(2))/2;(nd(1)+nd(2))/2;(wn(1)+wn(2))/2; (wp(1)+wp(2))/2 ];
%PP = @(p)-P(p(1),p(2),p(3),p(4),p(5));
%lb = [V(1),na(1),nd(1),wn(1),wp(1)];
%ub = [V(2),na(2),nd(2),wn(2),wp(2)];
-PP(P0)

ans = 
   0.012999986747818

%% Optimization
mx = fmincon(PP,P0,[],[],[],[],lb,ub);

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

-PP(mx)

ans = 
   0.024608689523057

V = mx(1)

V = 
   0.496682392530346

wn = mx(2)

wn = 
     4.000455217613459e-05

wp = mx(3)

wp = 
   0.015183822985277

nd = mx(4)

nd = 
     2.505000000000000e+20

na = mx(5)

na = 
     5.050000000000000e+15

%V = mx(1)
%na = mx(2)
%nd = mx(3)
%wn = mx(4)
%wp = mx(5)

[ANTONIO MANUEL]

You are a legend. Literally have no idea how much you just saved me and my partner. THANK YOU!!!