La meilleure manière de résoudre l'optimisation à plusieurs variables dans Matlab?

Je suis en train de calculer numériquement les solutions pour un système de plusieurs équations et des variables (+de 100). J'ai essayé jusqu'à présent, trois choses:

1. Je maintenant que le vecteur de p(i) (qui contient la plupart des variables endogènes) est décroissante. Donc j'ai donné simplement quelques points de départ, puis a été à la hausse(à la baisse) je suppose que quand j'ai vu que le p était trop faible(élevé). Bien sûr, cela a toujours été subordonnée à l'autre étant fixes, ce qui n'est pas le cas. Cela devrait finalement de travail, mais il n'est ni efficace, ni évident que j'en arrive à une solution en un temps fini. Il a travaillé lors de la réduction du système à 4 à 6 variables si.
2. Je pouvais créer 100+ boucles autour de l'autre et de l'utilisation non-bloquante pour chaque boucle. Ce serait éventuellement me conduire à la solution, mais de prendre des âges à la fois le programme (que je n'ai aucune idée de comment créer n boucles autour de l'autre sans avoir à écrire les boucles - qui est aussi mauvais que je le voudrais pour augmenter/diminuer la quantité de variables facilement) et à exécuter.
3. J'essayais fminsearch, mais comme prévu pour que étais quantité de variables - pas moyen!

Je vous serais reconnaissant de toutes les idées... Voici le code (celui de la fminsearch j'ai essayé):

C'est le fichier exécutable:

clear all
clc

% parameter

z=1.2;
w=20;
lam=0.7;
tau=1;
N=1000;
t_min=1;
t_max=4;
M=6;
a_min=0.6;
a_max=0.8;

t=zeros(1,N);
alp=zeros(1,M);
p=zeros(1,M);
p_min=2;
p_max=1;

for i=1:N
t(i)= t_min + (i-1)*(t_max - t_min)/(N-1);
end

for i=1:M
alp(i)= a_min + (i-1)*(a_max - a_min)/(M-1);
p(i)= p_min + (i-1)*(p_max - p_min)/(M-1);
end

fun=@(p) david(p ,z,w,lam,tau,N,M,t,alp);

p0=p;

fminsearch(fun,p0)

Et c'est le programme-fichier:

function crit=david(p, z,w,lam,tau,N,M,t,alp)
X = zeros(M,N);
pi = zeros(M,N);
C = zeros(1,N);
Xa=zeros(1,N);
Z=zeros(1,M);
rl=0.01;
rh=1.99;
EXD=140;
while (abs(EXD)>100)
r1=rl + 0.5*(rh-rl);  
for i=1:M
for j=1:N
X(i,j)=min(w*(1+lam), (alp(i) * p(i) / r1)^(1/(1-alp(i))) * t(j)^((z-alp(i))/(1-alp(i))));
pi(i,j)=p(i) * t(j)^(z-alp(i)) * X(i,j)^(alp(i)) - r1*X(i,j);
end
end
[C,I] = max(pi);
Xa(1)=X(I(1),1);
for j=2:N
Xa(j)=X(I(j),j);
end 
EXD=sum(Xa)- N*w;
if (abs(EXD)>100 && EXD>0)
rl=r1;
elseif (abs(EXD)>100 && EXD<0)
rh=r1;
end
end
Ya=zeros(M,N);
for j=1:N    
Ya(I(j),j)=t(j)^(z-alp(I(j))) * X(I(j),j)^(alp(I(j)));
end
Yi=sum(Ya,2);
if (Yi(1)==0)
Z(1)=-50;
end
for j=2:M
if (Yi(j)==0)
Z(j)=-50;
else
Z(j)=(p(1)/p(j))^tau - Yi(j)/Yi(1);
end
end
zz=sum(abs(Z))
crit=(sum(abs(Z)));