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