c0 = zeros(IL, JL, K, M);
p0 = ones(JW , K);
ptilda0 = zeros(IL, M, JL, K);
v0 = cell(IL, JL, JL, K, M);
iteration = 1;
epsilon = 1e-2;
while 1
c = sym('c', [IL, JL, K, M)];
objfun_c_symbExpression = f(c, p0, ptilda0, v0);
objfun_c_AnonymousFunction = matlabFunction(objfun_c_symbExpression, 'Vars', {c});
objfun_c = @(c) objfun_c_AnonymousFunction(c);
opts = optimset('Display','iter','Algorithm','sqp','MaxFunEval',inf,'MaxIter',Inf);
[c_opt,fval_c,flag_c] = fmincon(objfun_c,c0,A_c,b_c,[],[],lb_c,ub_c,[],opts);
p = sym('p', [JW , K]);
objfun_p_symbExpression = f(c_opt, p, ptilda0, v0);
objfun_p_AnonymousFunction = matlabFunction(objfun_p_symbExpression, 'Vars', {p});
objfun_p = @(p) objfun_p_AnonymousFunction(p);
opts = optimset('Display','iter','Algorithm','sqp','MaxFunEval',inf,'MaxIter',Inf);
[p_opt,fval_p,flag_p] = fmincon(objfun_p,p0,[],[],[],[],lb_p,ub_p,nonlcon_p,opts);
ptilda = sym('ptilda', [IL, M, JL, K]);
objfun_ptilda_symbExpression = f(c_opt, p_opt, ptilda, v0);
objfun_ptilda_AnonymousFunction = matlabFunction(objfun_ptilda_symbExpression, 'Vars', {ptilda});
objfun_ptilda = @(ptilda) objfun_ptilda_AnonymousFunction(ptilda);
opts = optimset('Display','iter','Algorithm','sqp','MaxFunEval',inf,'MaxIter',Inf);
[ptilda_opt,fval_ptilda,flag_ptilda] = fmincon(objfun_ptilda,ptilda0,A_ptilda,b_ptilda,[],[],lb_ptilda,ub_ptilda,nonlcon_ptilda,opts);
v = sym('v', ...
objfun_v_symbExpression = f(c_opt, p_opt, ptilda_opt, v);
objfun_v_AnonymousFunction = matlabFunction(objfun_v_symbExpression, 'Vars', {v});
objfun_v = @(v) objfun_v_AnonymousFunction(v);
opts = optimset('Display','iter','Algorithm','sqp','MaxFunEval',inf,'MaxIter',Inf);
[v_opt,fval_v,flag_v] = fmincon(objfun_v,v0,A_v,b_v,[],[],lb_v,ub_v,nonlcon_v,opts);
f_evolution = [f_evolution, f(c_opt, p_opt, ptilda_opt, v_opt)];
if size(f_evolution, 2)>2
err = abs(f_evolution(end) - f_evolution(end - 1) );
if err < epsilon
break,
end
end
iteration=iteration+1;
end