[x1,x2]=meshgrid(0:0.2:1,0:0.2:1); x3=1-x1-x2; mesh(x1,x2,x3) [n,m]=size(x3); for i=1:n*m; if x3(i)<0, x3(i)=0, end; end x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 -0.0000 -0.2000 0.6000 0.4000 0.2000 0 -0.2000 -0.4000 0.4000 0.2000 0 -0.2000 -0.4000 -0.6000 0.2000 0 -0.2000 -0.4000 -0.6000 -0.8000 0 0 -0.4000 -0.6000 -0.8000 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 -0.0000 -0.2000 0.6000 0.4000 0.2000 0 -0.2000 -0.4000 0.4000 0.2000 0 -0.2000 -0.4000 -0.6000 0.2000 0 0 -0.4000 -0.6000 -0.8000 0 0 -0.4000 -0.6000 -0.8000 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 -0.0000 -0.2000 0.6000 0.4000 0.2000 0 -0.2000 -0.4000 0.4000 0.2000 0 -0.2000 -0.4000 -0.6000 0.2000 0 0 -0.4000 -0.6000 -0.8000 0 0 0 -0.6000 -0.8000 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 -0.0000 -0.2000 0.6000 0.4000 0.2000 0 -0.2000 -0.4000 0.4000 0.2000 0 0 -0.4000 -0.6000 0.2000 0 0 -0.4000 -0.6000 -0.8000 0 0 0 -0.6000 -0.8000 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 -0.0000 -0.2000 0.6000 0.4000 0.2000 0 -0.2000 -0.4000 0.4000 0.2000 0 0 -0.4000 -0.6000 0.2000 0 0 0 -0.6000 -0.8000 0 0 0 -0.6000 -0.8000 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 -0.0000 -0.2000 0.6000 0.4000 0.2000 0 -0.2000 -0.4000 0.4000 0.2000 0 0 -0.4000 -0.6000 0.2000 0 0 0 -0.6000 -0.8000 0 0 0 0 -0.8000 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 0 -0.2000 0.6000 0.4000 0.2000 0 -0.2000 -0.4000 0.4000 0.2000 0 0 -0.4000 -0.6000 0.2000 0 0 0 -0.6000 -0.8000 0 0 0 0 -0.8000 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 0 -0.2000 0.6000 0.4000 0.2000 0 0 -0.4000 0.4000 0.2000 0 0 -0.4000 -0.6000 0.2000 0 0 0 -0.6000 -0.8000 0 0 0 0 -0.8000 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 0 -0.2000 0.6000 0.4000 0.2000 0 0 -0.4000 0.4000 0.2000 0 0 0 -0.6000 0.2000 0 0 0 -0.6000 -0.8000 0 0 0 0 -0.8000 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 0 -0.2000 0.6000 0.4000 0.2000 0 0 -0.4000 0.4000 0.2000 0 0 0 -0.6000 0.2000 0 0 0 0 -0.8000 0 0 0 0 -0.8000 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 0 -0.2000 0.6000 0.4000 0.2000 0 0 -0.4000 0.4000 0.2000 0 0 0 -0.6000 0.2000 0 0 0 0 -0.8000 0 0 0 0 0 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 0 0 0.6000 0.4000 0.2000 0 0 -0.4000 0.4000 0.2000 0 0 0 -0.6000 0.2000 0 0 0 0 -0.8000 0 0 0 0 0 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 0 0 0.6000 0.4000 0.2000 0 0 0 0.4000 0.2000 0 0 0 -0.6000 0.2000 0 0 0 0 -0.8000 0 0 0 0 0 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 0 0 0.6000 0.4000 0.2000 0 0 0 0.4000 0.2000 0 0 0 0 0.2000 0 0 0 0 -0.8000 0 0 0 0 0 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 0 0 0.6000 0.4000 0.2000 0 0 0 0.4000 0.2000 0 0 0 0 0.2000 0 0 0 0 0 0 0 0 0 0 -1.0000 x3 = 1.0000 0.8000 0.6000 0.4000 0.2000 0 0.8000 0.6000 0.4000 0.2000 0 0 0.6000 0.4000 0.2000 0 0 0 0.4000 0.2000 0 0 0 0 0.2000 0 0 0 0 0 0 0 0 0 0 0 mesh(x1,x2,x3) A=[1 1 1] A = 1 1 1 P=eye(3) - A'*inv(A*A')*A P = 0.6667 -0.3333 -0.3333 -0.3333 0.6667 -0.3333 -0.3333 -0.3333 0.6667 c=[1;2;3] c = 1 2 3 d=-c d = -1 -2 -3 dp=P*d dp = 1.0000 -0.0000 -1.0000 x0=[1/3 1/3 1/3]' x0 = 0.3333 0.3333 0.3333 a=-x0(3)/dp(3) a = 0.3333 gama=0.98; am=gama*a am = 0.3267 x1=x0+am*dp x1 = 0.6600 0.3333 0.0067 D=diag(x1) D = 0.6600 0 0 0 0.3333 0 0 0 0.0067 A1=A*D A1 = 0.6600 0.3333 0.0067 xk1=x1; [x1,x2]=meshgrid(0:0.2:1,0:0.2:1); x3=(1-A1(1)*x1-A1(2)*x2)/A1(3); for i=1:n*m; if x3(i)<0, x3(i)=0; end; end hold on mesh(x1,x2,x3) c1=(c'*D)' c1 = 0.6600 0.6667 0.0200 c c = 1 2 3 P=eye(3) - A1'*inv(A1*A1')*A1 P = 0.2033 -0.4024 -0.0080 -0.4024 0.7968 -0.0041 -0.0080 -0.0041 0.9999 dp1=P*(-c1) dp1 = 0.1342 -0.2655 -0.0120 -xk1(2)/dp1(2) ans = 1.2553 -xk1(3)/dp1(3) ans = 0.5566 xk1(3)+1.2553*dp1(3) ans = -0.0084 type karmarkar function [xo,Ot,nS,lambda]=karmarkar(c,x0,A,b,tol,mxit) % Solves the linear programming problem % min S(x)=c'x subject to: A x = b, x >= 0 % using the Karmarkar's algorithm. % % [xo,Ot,nS,lambda]=karmarkar(c,x0,A,b,tol,mxit) % % c: coefficients of the objective function % x0: initial point (must be feasible) % A: constraint matrix % b: constraint vector % tol: tolerance (default = 1e-4) % mxit: maximum number of iterations (default = 100) % xo: optimal point % Ot: optimal value of S % nS: number of objective function evaluations % lambda: Lagrange multipliers % Copyright (c) 2001 by LASIM-DEQUI-UFRGS % $Revision: 1.0 $ $Date: 2001/08/15 23:05:15 $ % Argimiro R. Secchi (arge@enq.ufrgs.br) if nargin < 4, error('karmarkar requires four input arguments'); end if nargin < 5 | isempty(tol), tol=1e-4; end if nargin < 6 | isempty(mxit), mxit=100; end [m,n]=size(A); erro=norm(c); I=eye(n); um=ones(n,1); gama=0.9995; nS=0; Q=I; while erro > tol & nS < mxit, nS=nS+1; B=inv(A*A')*A; P=I-A'*B; dp=-P*c; lambda=-B*c; z=-x0./dp; j=find(z>0); alfa=gama*min(z(j)); x1=x0+alfa*dp; xo=Q*x1 x01=Q*x0; %plot3([x01(1);xo(1)],[x01(2);xo(2)],[x01(3);xo(3)],'ro') D=diag(x1); Q=Q*D; A=A*D; c=D*c; x0=um; erro=norm(dp); end Ot=c'*xo; type ex_karma % Exemplo do algoritmo de Karmarkar % Edgar & Himmelblau, 1988 % % min c'x % A x = b % x >= 0 c=[1; 2; 3]; A=[1 1 1]; b=1; x0=[1/3; 1/3; 1/3] [xo,Ot,nS,lambda]=karmarkar(c,x0,A,b) c=[1; 2; 3]; A=[1 1 1]; b=1; x0=[1/3; 1/3; 1/3] x0 = 0.3333 0.3333 0.3333 [xo,Ot,nS,lambda]=karmarkar(c,x0,A,b) xo = 0.6665 0.3333 0.0002 xo = 0.9997 0.0002 0.0002 xo = 0.9999 0.0001 0.0000 xo = 1.0000 0.0000 0.0000 xo = 1.0000 0.0000 0.0000 Ot = 1.0000 nS = 4 lambda = -1.0000 edit karmarkar disp(sprintf('%d',feature('SessionTool'))) 0 ; dbstatus dbstack dbstack ; disp(which('karmarkar')); v:\cursos\pos\otimiza\aulas\karmarkar.m mdbstatus 'v:\cursos\pos\otimiza\aulas\karmarkar.m' eval( 'if exist(''must'',''var''), mauifunc(must), end', ''); eval( 'if exist(''mxit'',''var''), mauifunc(mxit), end', ''); eval( 'if exist(''size'',''var''), mauifunc(size), end', ''); disp(which('karmarkar')); v:\cursos\pos\otimiza\aulas\karmarkar.m clear 'v:\cursos\pos\otimiza\aulas\karmarkar.m' eval( 'if exist(''size'',''var''), mauifunc(size), end', ''); clear functions [xo,Ot,nS,lambda]=karmarkar(c,x0,A,b) xo = 0.6600 0.3333 0.0067 xo = 0.9870 0.0067 0.0064 xo = 0.9966 0.0032 0.0001 xo = 0.9998 0.0001 0.0001 xo = 1.0000 0.0000 0.0000 xo = 1.0000 0.0000 0.0000 xo = 1.0000 0.0000 0.0000 Ot = 1.0000 nS = 6 lambda = -1.0000 eval( 'if exist(''n'',''var''), mauifunc(n), end', ''); 6 eval( 'if exist(''norm'',''var''), mauifunc(norm), end', ''); eval( 'if exist(''n'',''var''), mauifunc(n), end', ''); 6 disp(which('karmarkar')); v:\cursos\pos\otimiza\aulas\karmarkar.m clear 'v:\cursos\pos\otimiza\aulas\karmarkar.m' clear functions Xp=eye(3) Xp = 1 0 0 0 1 0 0 0 1 x0=(Xp(1,:)+Xp(2,:)+Xp(3,:))/3 x0 = 0.3333 0.3333 0.3333 [xo,Ot,nS,lambda]=karmarkar(c,[0.5 0.5 0.5]',A,b) xo = 0.9997 0.5000 0.0002 xo = 1.4995 0.0003 0.0002 xo = 1.4999 0.0001 0.0000 xo = 1.5000 0.0000 0.0000 xo = 1.5000 0.0000 0.0000 xo = 1.5000 0.0000 0.0000 Ot = 2.2500 nS = 5 lambda = -1.0000 [xo,Ot,nS,lambda]=karmarkar(c,[0.5 0.5 0.5]'/norm([0.5 0.5 0.5]),A,b) xo = 1.1544 0.5774 0.0003 xo = 1.7315 0.0003 0.0003 xo = 1.7319 0.0001 0.0000 xo = 1.7321 0.0000 0.0000 xo = 1.7321 0.0000 0.0000 xo = 1.7321 0.0000 0.0000 Ot = 3.0000 nS = 5 lambda = -1.0000 x0=[0.5 0.5 0.5]' x0 = 0.5000 0.5000 0.5000 x0=x0/norm(x0) x0 = 0.5774 0.5774 0.5774 [xo,Ot,nS,lambda]=karmarkar(c,[0.5 0.5 0],A,b) ??? Error using ==> ./ Matrix dimensions must agree. Error in ==> v:\cursos\pos\otimiza\aulas\karmarkar.m On line 47 ==> z=-x0./dp; [xo,Ot,nS,lambda]=karmarkar(c,[0.5 0.4 0.1],A,b) ??? Error using ==> ./ Matrix dimensions must agree. Error in ==> v:\cursos\pos\otimiza\aulas\karmarkar.m On line 47 ==> z=-x0./dp; [xo,Ot,nS,lambda]=karmarkar(c,[0.5 0.4 0.1]',A,b) xo = 0.5999 0.4000 0.0000 xo = 0.9998 0.0002 0.0000 xo = 1.0000 0.0000 0.0000 xo = 1.0000 0.0000 0.0000 xo = 1.0000 0.0000 0.0000 Ot = 1.0000 nS = 4 lambda = -1.0000 type ex2_karm % Exemplo do algoritmo de Karmarkar % Edgar & Himmelblau, 1988 % % min c'x % A x = b % x >= 0 c0=[8.1; 10.8]; A0=[0.8 0.44 0.05 0.1 0.1 0.36]; b=[24000; 2000; 6000]; x0=[10000; 3000]; % incluindo as variaveis de folga: x0=[x0; b-A0*x0]; c=[-c0; zeros(3,1)]; A=[A0 eye(3)]; [xo,Ot,nS,lambda]=karmarkar(c,x0,A,b) xo=xo(1:2) Ot=c0'*xo ex2_karm xo = 1.0e+004 * 1.4327 1.2681 0.6958 0.0015 0.0002 xo = 1.0e+004 * 1.5025 1.2488 0.6486 0.0000 0.0002 xo = 1.0e+004 * 1.5970 1.2015 0.5937 0.0000 0.0078 xo = 1.0e+004 * 2.6202 0.6899 0.0003 0.0000 0.0896 xo = 1.0e+004 * 2.6207 0.6897 0.0000 0.0000 0.0897 xo = 1.0e+004 * 2.6207 0.6897 0.0000 0.0000 0.0897 xo = 1.0e+004 * 2.6207 0.6897 0.0000 0.0000 0.0897 xo = 1.0e+004 * 2.6191 0.6895 0.0000 0.0000 0.0897 xo = 1.0e+004 * 2.6191 0.6895 0.0000 0.0000 0.0897 Ot = -6.0698e+009 nS = 8 lambda = 4.6552 87.5172 0.0000 xo = 1.0e+004 * 2.6191 0.6895 Ot = 2.8661e+005 help linprog LINPROG Linear programming. X=LINPROG(f,A,b) solves the linear programming problem: min f'*x subject to: A*x <= b x X=LINPROG(f,A,b,Aeq,beq) solves the problem above while additionally satisfying the equality constraints Aeq*x = beq. X=LINPROG(f,A,b,Aeq,beq,LB,UB) defines a set of lower and upper bounds on the design variables, X, so that the solution is in the range LB <= X <= UB. Use empty matrices for LB and UB if no bounds exist. Set LB(i) = -Inf if X(i) is unbounded below; set UB(i) = Inf if X(i) is unbounded above. X=LINPROG(f,A,b,Aeq,beq,LB,UB,X0) sets the starting point to X0. This option is only available with the active-set algorithm. The default interior point algorithm will ignore any non-empty starting point. X=LINPROG(f,A,b,Aeq,Beq,LB,UB,X0,OPTIONS) minimizes with the default optimization parameters replaced by values in the structure OPTIONS, an argument created with the OPTIMSET function. See OPTIMSET for details. Use options are Display, Diagnostics, TolFun, LargeScale, MaxIter. Currently, only 'final' and 'off' are valid values for the parameter Display when LargeScale is 'off' ('iter' is valid when LargeScale is 'on'). [X,FVAL]=LINPROG(f,A,b) returns the value of the objective function at X: FVAL = f'*X. [X,FVAL,EXITFLAG] = LINPROG(f,A,b) returns EXITFLAG that describes the exit condition of LINPROG. If EXITFLAG is: > 0 then LINPROG converged with a solution X. 0 then LINPROG reached the maximum number of iterations without converging. < 0 then the problem was infeasible or LINPROG failed. [X,FVAL,EXITFLAG,OUTPUT] = LINPROG(f,A,b) returns a structure OUTPUT with the number of iterations taken in OUTPUT.iterations, the type of algorithm used in OUTPUT.algorithm, the number of conjugate gradient iterations (if used) in OUTPUT.cgiterations. [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=LINPROG(f,A,b) returns the set of Lagrangian multipliers LAMBDA, at the solution: LAMBDA.ineqlin for the linear inequalities A, LAMBDA.eqlin for the linear equalities Aeq, LAMBDA.lower for LB, and LAMBDA.upper for UB. NOTE: the LargeScale (the default) version of LINPROG uses a primal-dual method. Both the primal problem and the dual problem must be feasible for convergence. Infeasibility messages of either the primal or dual, or both, are given as appropriate. The primal problem in standard form is min f'*x such that A*x = b, x >= 0. The dual problem is max b'*y such that A'*y + s = f, s >= 0. op=optimset('LargeScale', 'off'); [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=(-c0,A0,b,[],[],[0 0],[],[],op) ??? AG,OUTPUT,LAMBDA]=( | Improper matrix assignment. A variable cannot return multiple arguments. [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=linprog(-c0,A0,b,[],[],[0 0],[],[],op) Optimization terminated successfully. X = 1.0e+004 * 2.6207 0.6897 FVAL = -2.8676e+005 EXITFLAG = 1 OUTPUT = iterations: 3 algorithm: 'medium-scale: activeset' firstorderopt: [] cgiterations: [] LAMBDA = lower: [2x1 double] upper: [2x1 double] eqlin: [0x1 double] ineqlin: [3x1 double] X X = 1.0e+004 * 2.6207 0.6897 xo xo = 1.0e+004 * 2.6191 0.6895 LAMBDA.ineqlin ans = 4.6552 87.5172 0 lambda lambda = 4.6552 87.5172 0.0000 op=optimset('LargeScale', 'on'); [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=linprog(-c0,A0,b,[],[],[0 0],[],x0,op) Warning: Interior Point method is ignoring starting point > In C:\Apps\Matlab\toolbox\optim\linprog.m at line 148 Optimization terminated successfully. X = 1.0e+004 * 2.6207 0.6897 FVAL = -2.8676e+005 EXITFLAG = 1 OUTPUT = iterations: 5 cgiterations: 0 algorithm: 'lipsol' LAMBDA = ineqlin: [3x1 double] eqlin: [0x1 double] upper: [2x1 double] lower: [2x1 double] what M-files in the current directory v:\cursos\pos\otimiza\aulas CATALIS ex_karma gtest13 nlconst test14 EXTRATOR ex_qp1 gtest2 nlp_internal test15 LUCRO ex_qp2 gtest9 powell test16 MINQUA ex_qp3 hkjeeves qpsub test17 MODELO ex_swarm htest1 refino test18 OPT_RES fmincon1 htest2 restr test19 PLANOS fminunc1 interior restr1 test2 READ2 fminusub karmarkar restr14 test20 SEMIDEF fun lmarqua restr15 test3 SMODELO gmilp1 lp_nlp restr16 test4 aurea gminlp1 milp1 restr17 test5 bandem1 gminlp2 minlp restr20 test6 bfgs gminlp3 minlp1 rosembr test7 bracket gminlp4 minlp2 set1 test8 buscarnd gminlp5 minlp3 setoptim test9 cgrad gminlp6 minlp4 sqp univar checkbounds gmodelagem minlp5 steepdes varmetr coggins gmurray minlp6 swarm writearq compdir grad minlpn test1 xplot complex grg modelagem test10 dfp gtest1 newton test11 dual gtest10 newton_h test12 ex2_karm gtest12 newtont test13 type ex_qp1 % min S(x) = 4 * x1^2 + 5 * x2^2 % s.a. h(x) = 2 * x1 + 3 * x2 - 6 = 0 % % x' = [1.071 1.286] Q=[4 0;0 5]*2 c=[0; 0] A=[] b=[] Aeq=[2 3] beq=6 lb=[] ub=[] op=optimset('LargeScale','off'); [x,s,ex,out,lambda]=quadprog(Q,c,A,b,Aeq,beq,lb,ub,[],op) ex_qp1 Q = 8 0 0 10 c = 0 0 A = [] b = [] Aeq = 2 3 beq = 6 lb = [] ub = [] Optimization terminated successfully. x = 1.0714 1.2857 s = 12.8571 ex = 1 out = iterations: 1 algorithm: 'medium-scale: active-set' firstorderopt: [] cgiterations: [] lambda = lower: [2x1 double] upper: [2x1 double] eqlin: -4.2857 ineqlin: [0x1 double] type ex_qp2 % min S(x) = (x1^2 + x2^2 + x3^2)/2 % s.a. g1(x) = 3/2 - x1 - x2 <= 0 % g2(x) = x1 - x3 -1 <= 0 % g3(x) = x2 + x3 - 2/5 <= 0 % % x' = [0.833 0.633 -0.2] (melhor solucao inviavel, no sentido de minimos quadrados) % x1' = [0.7667 0.5833 -0.1833] Q=[1 0 0;0 1 0;0 0 1] c=[0; 0; 0] A=[-1 -1 0 1 0 -1 0 1 1] b=[-3/2; 1; 2/5] Aeq=[] beq=[] lb=[] ub=[] [x,s,ex,out,lambda]=quadprog(Q,c,A,b,Aeq,beq,lb,ub) b1=[-3/2*0.9; 1; 2/5] [x1,s1,ex,out,lambda1]=quadprog(Q,c,A,b1,Aeq,beq,lb,ub) Q=[1 0 0;0 1 0;0 0 1] Q = 1 0 0 0 1 0 0 0 1 c=[0; 0; 0] c = 0 0 0 A=[-1 -1 0 1 0 -1 0 1 1] A = -1 -1 0 1 0 -1 0 1 1 b=[-3/2; 1; 2/5] b = -1.5000 1.0000 0.4000 Aeq=[] Aeq = [] beq=[] beq = [] lb=[] lb = [] ub=[] ub = [] [x,s,ex,out,lambda]=quadprog(Q,c,A,b,Aeq,beq,lb,ub) Exiting: The constraints are overly stringent; no feasible starting point found. x = 0.8333 0.6333 -0.2000 s = 0.5678 ex = -1 out = iterations: 0 algorithm: 'medium-scale: active-set' firstorderopt: [] cgiterations: [] lambda = lower: [3x1 double] upper: [3x1 double] eqlin: [0x1 double] ineqlin: [3x1 double] A A = -1 -1 0 1 0 -1 0 1 1 inv(A)*b Warning: Matrix is singular to working precision. ans = NaN NaN NaN type ex_qp2 % min S(x) = (x1^2 + x2^2 + x3^2)/2 % s.a. g1(x) = 3/2 - x1 - x2 <= 0 % g2(x) = x1 - x3 -1 <= 0 % g3(x) = x2 + x3 - 2/5 <= 0 % % x' = [0.833 0.633 -0.2] (melhor solucao inviavel, no sentido de minimos quadrados) % x1' = [0.7667 0.5833 -0.1833] Q=[1 0 0;0 1 0;0 0 1] c=[0; 0; 0] A=[-1 -1 0 1 0 -1 0 1 1] b=[-3/2; 1; 2/5] Aeq=[] beq=[] lb=[] ub=[] [x,s,ex,out,lambda]=quadprog(Q,c,A,b,Aeq,beq,lb,ub) b1=[-3/2*0.9; 1; 2/5] [x1,s1,ex,out,lambda1]=quadprog(Q,c,A,b1,Aeq,beq,lb,ub) LAMBDA LAMBDA = ineqlin: [3x1 double] eqlin: [0x1 double] upper: [2x1 double] lower: [2x1 double] LAMBDA.ineqlin ans = 4.6552 87.5172 0.0000 type ex_qp2 % min S(x) = (x1^2 + x2^2 + x3^2)/2 % s.a. g1(x) = 3/2 - x1 - x2 <= 0 % g2(x) = x1 - x3 -1 <= 0 % g3(x) = x2 + x3 - 2/5 <= 0 % % x' = [0.833 0.633 -0.2] (melhor solucao inviavel, no sentido de minimos quadrados) % x1' = [0.7667 0.5833 -0.1833] Q=[1 0 0;0 1 0;0 0 1] c=[0; 0; 0] A=[-1 -1 0 1 0 -1 0 1 1] b=[-3/2; 1; 2/5] Aeq=[] beq=[] lb=[] ub=[] [x,s,ex,out,lambda]=quadprog(Q,c,A,b,Aeq,beq,lb,ub) b1=[-3/2*0.9; 1; 2/5] [x1,s1,ex,out,lambda1]=quadprog(Q,c,A,b1,Aeq,beq,lb,ub) lamdba ??? Undefined function or variable 'lamdba'. lambda lambda = lower: [3x1 double] upper: [3x1 double] eqlin: [0x1 double] ineqlin: [3x1 double] lambda.ineqlin ans = 0.2357 0.2357 0.2357 A*x ans = -1.4667 1.0333 0.4333 b b = -1.5000 1.0000 0.4000 type ex_qp2 % min S(x) = (x1^2 + x2^2 + x3^2)/2 % s.a. g1(x) = 3/2 - x1 - x2 <= 0 % g2(x) = x1 - x3 -1 <= 0 % g3(x) = x2 + x3 - 2/5 <= 0 % % x' = [0.833 0.633 -0.2] (melhor solucao inviavel, no sentido de minimos quadrados) % x1' = [0.7667 0.5833 -0.1833] Q=[1 0 0;0 1 0;0 0 1] c=[0; 0; 0] A=[-1 -1 0 1 0 -1 0 1 1] b=[-3/2; 1; 2/5] Aeq=[] beq=[] lb=[] ub=[] [x,s,ex,out,lambda]=quadprog(Q,c,A,b,Aeq,beq,lb,ub) b1=[-3/2*0.9; 1; 2/5] [x1,s1,ex,out,lambda1]=quadprog(Q,c,A,b1,Aeq,beq,lb,ub) b1=[-3/2*0.9; 1; 2/5] b1 = -1.3500 1.0000 0.4000 [x1,s1,ex,out,lambda1]=quadprog(Q,c,A,b1,Aeq,beq,lb,ub) Optimization terminated successfully. x1 = 0.7667 0.5833 -0.1833 s1 = 0.4808 ex = 1 out = iterations: 3 algorithm: 'medium-scale: active-set' firstorderopt: [] cgiterations: [] lambda1 = lower: [3x1 double] upper: [3x1 double] eqlin: [0x1 double] ineqlin: [3x1 double] A*x1 ans = -1.3500 0.9500 0.4000 b1 b1 = -1.3500 1.0000 0.4000 quit type ex_qp1 % min S(x) = 4 * x1^2 + 5 * x2^2 % s.a. h(x) = 2 * x1 + 3 * x2 - 6 = 0 % % x' = [1.071 1.286] Q=[4 0;0 5]*2 c=[0; 0] A=[] b=[] Aeq=[2 3] beq=6 lb=[] ub=[] op=optimset('LargeScale','off'); [x,s,ex,out,lambda]=quadprog(Q,c,A,b,Aeq,beq,lb,ub,[],op) type ex_qp2 % min S(x) = (x1^2 + x2^2 + x3^2)/2 % s.a. g1(x) = 3/2 - x1 - x2 <= 0 % g2(x) = x1 - x3 -1 <= 0 % g3(x) = x2 + x3 - 2/5 <= 0 % % x' = [0.833 0.633 -0.2] (melhor solucao inviavel, no sentido de minimos quadrados) % x1' = [0.7667 0.5833 -0.1833] Q=[1 0 0;0 1 0;0 0 1] c=[0; 0; 0] A=[-1 -1 0 1 0 -1 0 1 1] b=[-3/2; 1; 2/5] Aeq=[] beq=[] lb=[] ub=[] [x,s,ex,out,lambda]=quadprog(Q,c,A,b,Aeq,beq,lb,ub) b1=[-3/2*0.9; 1; 2/5] [x1,s1,ex,out,lambda1]=quadprog(Q,c,A,b1,Aeq,beq,lb,ub) type ex_qp1 % min S(x) = 4 * x1^2 + 5 * x2^2 % s.a. h(x) = 2 * x1 + 3 * x2 - 6 = 0 % % x' = [1.071 1.286] Q=[4 0;0 5]*2 c=[0; 0] A=[] b=[] Aeq=[2 3] beq=6 lb=[] ub=[] op=optimset('LargeScale','off'); [x,s,ex,out,lambda]=quadprog(Q,c,A,b,Aeq,beq,lb,ub,[],op) eig([8 2;2 10]0 ??? eig([8 2;2 10]0 | Improper function reference. A "," or ")" is expected. eig([8 2;2 10]0 ??? eig([8 2;2 10]0 | Improper function reference. A "," or ")" is expected. eig([8 2;2 10]) ans = 6.7639 11.2361 type ex_qp2 % min S(x) = (x1^2 + x2^2 + x3^2)/2 % s.a. g1(x) = 3/2 - x1 - x2 <= 0 % g2(x) = x1 - x3 -1 <= 0 % g3(x) = x2 + x3 - 2/5 <= 0 % % x' = [0.833 0.633 -0.2] (melhor solucao inviavel, no sentido de minimos quadrados) % x1' = [0.7667 0.5833 -0.1833] Q=[1 0 0;0 1 0;0 0 1] c=[0; 0; 0] A=[-1 -1 0 1 0 -1 0 1 1] b=[-3/2; 1; 2/5] Aeq=[] beq=[] lb=[] ub=[] [x,s,ex,out,lambda]=quadprog(Q,c,A,b,Aeq,beq,lb,ub) b1=[-3/2*0.9; 1; 2/5] [x1,s1,ex,out,lambda1]=quadprog(Q,c,A,b1,Aeq,beq,lb,ub) [x1,x2]=meshgrid(0:10,0:10); S=-6*x1-2*x2+2*x1.*x2+4*x1.^2+5*x2.^2; mesh(x1,x2,S) [x1,x2]=meshgrid(0:0.1:1,0:0.1:1); S=-6*x1-2*x2+2*x1.*x2+4*x1.^2+5*x2.^2; mesh(x1,x2,S) contour(x1,x2,S) hold Current plot held X1=0:0.1:1; X21=(8-2*X1)/8; X22=(10-4*X1); plot(X1,X21,'r',X1,X22,'b') [x1,x2]=meshgrid(0:10,0:10); S=-6*x1-2*x2+2*x1.*x2+4*x1.^2+5*x2.^2; X1=0:1:10; X21=(8-2*X1)/8; X22=(10-4*X1); contour(x1,x2,S) contour(x1,x2,S,50) hold Current plot held plot(X1,X21,'r',X1,X22,'b') eig([8 2;2 10]) ans = 6.7639 11.2361 plot(2.4,0.4,'*'); Q=[8 2;2 10]; c=[-6;-2]; xs=-Q*c xs = 52 32 xs=-inv(Q)*c xs = 0.7368 0.0526 plot(xs(1),xs(2),'*r') S=-20*x1-8*x2+2*x1.*x2+4*x1.^2+5*x2.^2; contour(x1,x2,S,50) S=-80*x1-40*x2+2*x1.*x2+4*x1.^2+5*x2.^2; contour(x1,x2,S,50) contour(x1,x2,S,50) hold Current plot held plot(X1,X21,'r',X1,X22,'b') c=[-80;-40] c = -80 -40 xs=-inv(Q)*c xs = 9.4737 2.1053 plot(xs(1),xs(2),'*r') plot(2.4,0.4,'*'); b=[8;10] b = 8 10 A=[2 8;4 1] A = 2 8 4 1 ch=0.5*[c' 0 0 -b' 0 0]' ch = -40 -20 0 0 -4 -5 0 0 bh=[b' -c']' bh = 8 10 80 40 Ah=[A eye(2) zeros(2) zeros(2);Q zeros(2) A' -eye(2)] Ah = 2 8 1 0 0 0 0 0 4 1 0 1 0 0 0 0 8 2 0 0 2 4 -1 0 2 10 0 0 8 1 0 -1 help linprog LINPROG Linear programming. X=LINPROG(f,A,b) solves the linear programming problem: min f'*x subject to: A*x <= b x X=LINPROG(f,A,b,Aeq,beq) solves the problem above while additionally satisfying the equality constraints Aeq*x = beq. X=LINPROG(f,A,b,Aeq,beq,LB,UB) defines a set of lower and upper bounds on the design variables, X, so that the solution is in the range LB <= X <= UB. Use empty matrices for LB and UB if no bounds exist. Set LB(i) = -Inf if X(i) is unbounded below; set UB(i) = Inf if X(i) is unbounded above. X=LINPROG(f,A,b,Aeq,beq,LB,UB,X0) sets the starting point to X0. This option is only available with the active-set algorithm. The default interior point algorithm will ignore any non-empty starting point. X=LINPROG(f,A,b,Aeq,Beq,LB,UB,X0,OPTIONS) minimizes with the default optimization parameters replaced by values in the structure OPTIONS, an argument created with the OPTIMSET function. See OPTIMSET for details. Use options are Display, Diagnostics, TolFun, LargeScale, MaxIter. Currently, only 'final' and 'off' are valid values for the parameter Display when LargeScale is 'off' ('iter' is valid when LargeScale is 'on'). [X,FVAL]=LINPROG(f,A,b) returns the value of the objective function at X: FVAL = f'*X. [X,FVAL,EXITFLAG] = LINPROG(f,A,b) returns EXITFLAG that describes the exit condition of LINPROG. If EXITFLAG is: > 0 then LINPROG converged with a solution X. 0 then LINPROG reached the maximum number of iterations without converging. < 0 then the problem was infeasible or LINPROG failed. [X,FVAL,EXITFLAG,OUTPUT] = LINPROG(f,A,b) returns a structure OUTPUT with the number of iterations taken in OUTPUT.iterations, the type of algorithm used in OUTPUT.algorithm, the number of conjugate gradient iterations (if used) in OUTPUT.cgiterations. [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=LINPROG(f,A,b) returns the set of Lagrangian multipliers LAMBDA, at the solution: LAMBDA.ineqlin for the linear inequalities A, LAMBDA.eqlin for the linear equalities Aeq, LAMBDA.lower for LB, and LAMBDA.upper for UB. NOTE: the LargeScale (the default) version of LINPROG uses a primal-dual method. Both the primal problem and the dual problem must be feasible for convergence. Infeasibility messages of either the primal or dual, or both, are given as appropriate. The primal problem in standard form is min f'*x such that A*x = b, x >= 0. The dual problem is max b'*y such that A'*y + s = f, s >= 0. op=optimset('LargeScale','off'); [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=LINPROG(ch,Ah,bh,[],[],zeros(8,1),[],[],op) Exiting: The solution is unbounded and at infinity; the constraints are not restrictive enough. X = 1.0e+015 * 0.0000 0.0000 -0.0000 0 0 2.3570 9.4281 2.3570 FVAL = -1.1785e+016 EXITFLAG = -1 OUTPUT = iterations: 8 algorithm: 'medium-scale: activeset' firstorderopt: [] cgiterations: [] LAMBDA = lower: [8x1 double] upper: [8x1 double] eqlin: [0x1 double] ineqlin: [4x1 double] [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=LINPROG(ch,[],[],Ah,bh,zeros(8,1),[],[],op) Exiting: The solution is unbounded and at infinity; the constraints are not restrictive enough. X = 1.0e+015 * 0.0000 0.0000 0.0000 -0.0000 0 2.3570 9.4281 2.3570 FVAL = -1.1785e+016 EXITFLAG = -1 OUTPUT = iterations: 4 algorithm: 'medium-scale: activeset' firstorderopt: [] cgiterations: [] LAMBDA = lower: [8x1 double] upper: [8x1 double] eqlin: [4x1 double] ineqlin: [0x1 double] c c = -80 -40 ch ch = -40 -20 0 0 -4 -5 0 0 bh bh = 8 10 80 40 A A = 2 8 4 1 Ah Ah = 2 8 1 0 0 0 0 0 4 1 0 1 0 0 0 0 8 2 0 0 2 4 -1 0 2 10 0 0 8 1 0 -1 zeros(8,1) ans = 0 0 0 0 0 0 0 0 [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=LINPROG(ch,[],[],Ah,bh,zeros(8,1),[],[],op) Exiting: The solution is unbounded and at infinity; the constraints are not restrictive enough. X = 1.0e+015 * 0.0000 0.0000 0.0000 -0.0000 0 2.3570 9.4281 2.3570 FVAL = -1.1785e+016 EXITFLAG = -1 OUTPUT = iterations: 4 algorithm: 'medium-scale: activeset' firstorderopt: [] cgiterations: [] LAMBDA = lower: [8x1 double] upper: [8x1 double] eqlin: [4x1 double] ineqlin: [0x1 double] LAMBDA.eqlin ans = 0 0 0 0 LAMBDA.upper ans = 0 0 0 0 0 0 0 0 LAMBDA.lower ans = 0 0 0 0 0 0 0 0 ch=0.5*[c' 0 0 b' 0 0]' ch = -40 -20 0 0 4 5 0 0 [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=LINPROG(ch,[],[],Ah,bh,zeros(8,1),[],[],op) Optimization terminated successfully. X = 2.4000 0.4000 0 0 2.1600 13.9200 0 0 FVAL = -25.7600 EXITFLAG = 1 OUTPUT = iterations: 4 algorithm: 'medium-scale: activeset' firstorderopt: [] cgiterations: [] LAMBDA = lower: [8x1 double] upper: [8x1 double] eqlin: [4x1 double] ineqlin: [0x1 double] [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=LINPROG(ch,[],[],Ah,bh,zeros(8,1),[],[],op) Optimization terminated successfully. X = 2.4000 0.4000 0 0 2.1600 13.9200 0 0 FVAL = -25.7600 EXITFLAG = 1 OUTPUT = iterations: 4 algorithm: 'medium-scale: activeset' firstorderopt: [] cgiterations: [] LAMBDA = lower: [8x1 double] upper: [8x1 double] eqlin: [4x1 double] ineqlin: [0x1 double] eye([0 1;1 0]) ??? Error using ==> eye Unknown command option. eig([0 1;1 0]) ans = -1.0000 1.0000 eig([1 1;1 1]) ans = 0 2.0000 eig([-1 1;1 -1]) ans = -2.0000 -0.0000