ObjFun = f( R0(c), R1(c));
Const1 = sum(sum(c(:,:,:,:))) - 1 <= 0;
for k = 1 : K
for m = 1 : M
for l = 1 : L
for i = 1 : I
Const2 = -1*c(i, l, k, m) <= 0;
Const3 = c(i, l, k, m) - 1 <= 0;
end
end
end
end
LagrangianFun = ObjFun + lambda1*Const1 + lambda2*Const2 + lambda3*Const3;
dLagrangianFun_dc = diff(LagrangianFun,l) == 0;
dLagrangianFun_dlambda1 = diff(LagrangianFun,lambda1) == 0;
dLagrangianFun_dlambda2 = diff(LagrangianFun,lambda2) == 0;
dLagrangianFun_dlambda3 = diff(LagrangianFun,lambda3) == 0;
system = [dLagrangianFun_dc; dLagrangianFun_dlambda1; dLagrangianFun_dlambda2; dLagrangianFun_dlambda3];
[c_val, lambda1_val, , lambda2_val, lambda3_val] = solve(system, [c lambda1 lambda2 lambda3 ], 'Real', true) ;