%Problema de optimizacion clear; clc; close all; syms x1 x2 u1 u2 u3 u4; f =(2*sin(1/6*x1)+5).*x1 + (2*sin(1/6*x2)+10).*x2 + 2*log10(x2) + log10(x1); g1 = -x1-x2+7.2E3; g2 = -x1+4*x2; g3 = -x1; g4 = -x2; %Lagrange L = f;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) sol1 = comprobacion(sol,g1,g2,g3,g4); disp(sol1); end L = f+u1*g1;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u1 >= 0 sol2 = comprobacion(sol,g1,g2,g3,g4); disp(sol2); end end L = f+u2*g2;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u2 >= 0 sol3 = comprobacion(sol,g1,g2,g3,g4); disp(sol2); end end L = f+u1*g1+u2*g2;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u1 >= 0 if sol.u2 >= 0 sol4 = comprobacion(sol,g1,g2,g3,g4); disp(sol4); end end end L = f+u3*g3;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u3 >= 0 sol5 = comprobacion(sol,g1,g2,g3,g4); disp(sol5) end end L = f+u1*g1+u3*g3;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u1 >= 0 if sol.u3 >= 0 sol6 = comprobacion(sol,g1,g2,g3,g4); disp(sol6); end end end L = f+u2*g2+u3*g3;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u2 >= 0 if sol.u3 >= 0 sol7 = comprobacion(sol,g1,g2,g3,g4); disp(sol7); end end end L = f+u1*g1+u2*g2+u3*g3;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u1 >= 0 if sol.u2 >= 0 if sol.u3 >= 0 sol8 = comprobacion(sol,g1,g2,g3,g4); disp(sol8); end end end end L = f+u4*g4;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u4 >= 0 sol9 = comprobacion(sol,g1,g2,g3,g4); disp(sol9); end end L = f+u1*g1+u4*g4;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u1 >= 0 if sol.u4 >= 0 sol10 = comprobacion(sol,g1,g2,g3,g4); disp(sol10); end end end L = f+u2*g2+u4*g4;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u2 >= 0 if sol.u4 >= 0 sol11 = comprobacion(sol,g1,g2,g3,g4); disp(sol11); end end end L = f+u1*g1+u2*g2+u4*g4;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u1 >= 0 if sol.u2 if sol.u4 >= 0 sol12 = comprobacion(sol,g1,g2,g3,g4); disp(sol12); end end end end L = f+u3*g3+u4*g4;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u3 >= 0 if sol.u4 >= 0 sol13 = comprobacion(sol,g1,g2,g3,g4); disp(sol13); end end end L = f+u1*g1+u3*g3+u4*g4;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u1 >= 0 if sol.u3 >= 0 if sol.u4 >= 0 sol14 = comprobacion(sol,g1,g2,g3,g4); disp(sol14); end end end end L = f+u2*g2+u3*g3+u4*g4;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u2 >= 0 if sol.u3 >= 0 if sol.u4 >= 0 sol15 = comprobacion(sol,g1,g2,g3,g4); disp(sol15); end end end end L = f+u1*g1+u2*g2+u3*g3+u4*g4;sol = solve(gradient(L)); if ~any(structfun(@isempty, sol)) if sol.u1 >= 0 if sol.u2 >= 0 if sol.u3 >= 0 if sol.u4 >= 0 sol16 = comprobacion(sol,g1,g2,g3,g4); disp(sol16); end end end end end %Dual----------------------------------------------- Lag = f+u1*g1+u2*g2+u3*g3+u4*g4; sDual = solve(gradient(Lag,[x1,x2])); q = subs(Lag,{x1 x2},{sDual.x1 sDual.x2}); if ~any(structfun(@isempty, sol)) s2Dual = solve(gradient(q)); if s2Dual.u1 >= 0 if s2Dual.u2 >= 0 if s2Dual.u3 >= 0 if s2Dual.u4 >= 0 sDual.x1=subs(sDual.x1,u1,s2Dual.u1(1));sDual.x2=subs(sDual.x2,u1,s2Dual.u1(1)); sDual.x1=subs(sDual.x1,u2,s2Dual.u2(1));sDual.x2=subs(sDual.x2,u2,s2Dual.u2(1)); sDual.x1=subs(sDual.x1,u3,s2Dual.u3(1));sDual.x2=subs(sDual.x2,u3,s2Dual.u3(1)); sDual.x1=subs(sDual.x1,u4,s2Dual.u4(1));sDual.x2=subs(sDual.x2,u4,s2Dual.u4(1)); disp(sDual); end end end end else disp("No se encuentra solución Dual") end %Penalizacion----------------------------------------------------- xk1 = [9000;500]; ck = 100; for i=1:0.5:50 %Fijo el valor de ck fun=@(x)funM1Pen(x,ck); %Optimizo con el punto y c actual xknew=fminunc(fun,xk1); %Compruebo si llego al fin [m,ff,g]=fun(xknew); xk1=xknew; if norm(g)<1e-5 break; end %Incremento ck ck=ck+1; xk1=xknew; %disp([ck,xk',g]) %Muestro end if subs(g1,{x1 x2},{xk1(1) xk1(2)}) <= 0 if subs(g2,{x1 x2},{xk1(1) xk1(2)}) <= 0 if subs(g3,{x1 x2},{xk1(1) xk1(2)}) <= 0 if subs(g4,{x1 x2},{xk1(1) xk1(2)}) <= 0 sprintf('Penalización -> x1 = %d x2 = %d',xk1') eval(subs(f,{x1 x2},{xk1(1) xk1(2)})) end end end end %Barrera --------------------------------------------- xk=[8000;400]; %PI(dentro de la zona) ck=0.1; %c_0: muy pequeÒo for i=1:40 %Fijo el valor de ck fun=@(x)funM1Bar(x,ck); %Optimizo con el punto y c actual xknew=fminsearch(fun,xk); %Compruebo si llego al fin [m,fb,g]=fun(xknew); %Incremento ck ck=ck*2; xk=xknew; %disp([ck,xk',m,f,g]) %Muestro if max(g) <= 0 break; end end if subs(g1,{x1 x2},{xk(1) xk(2)}) <= 0 if subs(g2,{x1 x2},{xk(1) xk(2)}) <= 0 if subs(g3,{x1 x2},{xk(1) xk(2)}) <= 0 if subs(g4,{x1 x2},{xk(1) xk(2)}) <= 0 sprintf("Barrera -> x1 = %d x2 = %d",xk') eval(subs(f,{x1 x2},{xk(1) xk(2)})) end end end end %Fmincon------------------------------------------ x0=[0,0]'; fun=@(x)funM1Eval(x); nonlcon=[]; opt=optimset('Display','iter','MaxIter',50); A=[-1,-1;-1 4];b=[-7.2E3;0]; %Lineales A*x≤b Aeq=[];beq=[];lb=[0;0];ub=[]; xk2 = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,opt); disp(xk2); eval(subs(f,{x1 x2},{xk2(1) xk2(2)})) grid_points = 100; maxX1 = 30000; maxX2 = 30000; xg = 0:0.1:30000; x1g = linspace(0,maxX1,grid_points); x2g = linspace(0,maxX2,grid_points); [x1g, x2g] = meshgrid(x1g,x2g); fg = (2*sin(1/6*x1g)+5).*x1g + (2*sin(1/6*x2g)+10).*x2g + 2*log10(x2g) + log10(x1g); g1g = 7.2E3-xg; g2g = xg./4; %g3g = -x1; %g4g = -x2; close all; figure(); contour(x1g,x2g,fg,ShowText="on"); hold on; plot(g1g,xg); plot(g2g,xg); plot(xk(1),xk(2),'*', 'markersize', 8); plot(xk1(1),xk1(2),'*', 'markersize', 8); plot(xk2(1),xk2(2),'*', 'markersize', 8); plot(sol4.x1,sol4.x2,'*', 'markersize', 8); xlim([0 30000]); function f=funM1Eval(x) x2=x(2); x1=x(1); f =(0.5*sin(1/6*x1)+5).*x1 + (5*sin(1/6*x2)+10).*x2 + 2*log10(x2) + log10(x1); %f = -f; end function [m,f,g]=funM1Bar(x,c) x2=x(2); x1=x(1); f = (2*sin(1/6*x1)+5).*x1 + (2*sin(1/6*x2)+10).*x2 + 2*log10(x2) + log10(x1); g(1)=-x1-x2+7.2E3; %g1(x)≤0 g(2)=-x1+4*x2; %g2(x)≤0 g(3)=-x1; g(4)=-x2; m=f-sum(log10(-g)); %m(x,c) end function [m,f,g] = funM1Pen(x,c) x2 = x(2); x1 = x(1); f =(2*sin(1/6*x1)+5).*x1 + (2*sin(1/6*x2)+10).*x2 + 2*log10(x2) + log10(x1); g(1) = max([0,-x1-x2+7.2E3]); g(2) = max([0,-x1+4*x2]); g(3) = max([0,-x1]); g(4) = max([0,-x2]); m=f+c*norm(g); end %