clc cd V:\Cursos\Pos\Otimiza\Aulas help grg Constrained optimization using GRG to solve problems like: min S(x) subject to: g(x) <= 0, h(x) = 0 (nonlinear constraints) x Lb <= x <= Ub [xo,Ot,nS,lambda]=grg(S,Res,x0,iXI,ip,GrS,GrR,Lb,Ub,problem,almax,opt) S: objective function x0: initial point ixI: index vector of independent variables (default: slack variables + first ones) ip: (0) no plot (default), (>0) plot figure ip with pause, (<0) plot figure ip Res: contraint function returning [g(x),h(x)] GrS: gradient vector function of S(x) GrR: gradient matrix function of Res(x), returning [Dg(x),Dh(x)] Lb, Ub: lower and upper bound vectors [used also to plot (plot default = x0*(1+/-2))] problem: (-1): minimum (default), (1): maximum almax: maximum stepsize (default = 10) opt: options vector set using optimset('grg') xo: optimal point Ot: optimal value of S nS: number of objective function evaluations lambda: Lagrange multipliers opt=optimset('grg') opt = ActiveConstrTol: [] DerivativeCheck: 'off' Diagnostics: 'off' DiffMaxChange: [] DiffMinChange: [] Display: 'final' GoalsExactAchieve: [] GradConstr: 'on' GradObj: 'on' Hessian: [] HessMult: [] HessPattern: [] HessUpdate: [] Jacobian: [] JacobMult: [] JacobPattern: [] LargeScale: [] LevenbergMarquardt: [] LineSearchType: [] MaxFunEvals: [] MaxIter: '50*(1+4*~(ip>0))' MaxPCGIter: [] MeritFunction: [] MinAbsMax: [] Preconditioner: [] PrecondBandWidth: [] ShowStatusWindow: [] TolCon: 1.0000e-006 TolFun: 1.0000e-006 TolPCG: [] TolX: 1.0000e-006 TypicalX: [] opt=optimset(opt,'MaxIter',200) opt = ActiveConstrTol: [] DerivativeCheck: 'off' Diagnostics: 'off' DiffMaxChange: [] DiffMinChange: [] Display: 'final' GoalsExactAchieve: [] GradConstr: 'on' GradObj: 'on' Hessian: [] HessMult: [] HessPattern: [] HessUpdate: [] Jacobian: [] JacobMult: [] JacobPattern: [] LargeScale: [] LevenbergMarquardt: [] LineSearchType: [] MaxFunEvals: [] MaxIter: 200 MaxPCGIter: [] MeritFunction: [] MinAbsMax: [] Preconditioner: [] PrecondBandWidth: [] ShowStatusWindow: [] TolCon: 1.0000e-006 TolFun: 1.0000e-006 TolPCG: [] TolX: 1.0000e-006 TypicalX: [] help grgh grgh.m not found. help grg Constrained optimization using GRG to solve problems like: min S(x) subject to: g(x) <= 0, h(x) = 0 (nonlinear constraints) x Lb <= x <= Ub [xo,Ot,nS,lambda]=grg(S,Res,x0,iXI,ip,GrS,GrR,Lb,Ub,problem,almax,opt) S: objective function x0: initial point ixI: index vector of independent variables (default: slack variables + first ones) ip: (0) no plot (default), (>0) plot figure ip with pause, (<0) plot figure ip Res: contraint function returning [g(x),h(x)] GrS: gradient vector function of S(x) GrR: gradient matrix function of Res(x), returning [Dg(x),Dh(x)] Lb, Ub: lower and upper bound vectors [used also to plot (plot default = x0*(1+/-2))] problem: (-1): minimum (default), (1): maximum almax: maximum stepsize (default = 10) opt: options vector set using optimset('grg') xo: optimal point Ot: optimal value of S nS: number of objective function evaluations lambda: Lagrange multipliers type test14 function S=test(x) % Edgar & Himmelblau, 1988 % x0 = [-1, -2]' % sem restricao: % xo = [2, 1]' % S(xo) = 0 % com restricao: % xo = [0.8229, 0.9114]' % S(xo) = 1.3935 S=(x(1)-2).^2+(x(2)-1).^2; type restr14 function [g,h]=restr(x) % Edgar & Himmelblau, 1988 % x0 = [-1, -2]' % xo = [0.8229, 0.9114]' % S(xo) = 1.3935 g=(x(1).^2)/4+x(2).^2-1; h=x(1)-2*x(2)+1; [xo,Ot,nS,lambda]=grg('test14','restr14',[-1 -2],1,1,'gtest14','grestr14',[-5 -5],[5 5]) Pause: hit any key to continue... max Directional Call Iter F-count f(x) constraint Step-size derivative 1 52 1 3.93 0.2 -81 0 Pause: hit any key to continue... 2 103 0.2 1.77 0.2 -4 0 Pause: hit any key to continue... 3 154 0.2 1.77 1.08e+006 -4 1 Pause: hit any key to continue... 4 204 0.2 1.77 5.39e+006 -4.57e-030 1 Pause: hit any key to continue... 5 282 0.2 1.77 1.32 -1.19e-016 1 Pause: hit any key to continue... 6 354 0.2 1.77 0.67 -5.97e-016 1 Pause: hit any key to continue... 7 414 0.2 1.77 1.04 -1.59e-015 1 Pause: hit any key to continue... 8 483 0.2 1.77 1.05 -6.29e-015 1 Pause: hit any key to continue... 9 546 0.2 1.77 1.06 -2.5e-014 1 Pause: hit any key to continue... 10 611 0.2 1.77 1.07 -1e-013 1 Pause: hit any key to continue... 11 682 0.2 1.77 1.06 -4.06e-013 1 Pause: hit any key to continue... 12 735 0.2 1.77 1.06 -1.62e-012 1 Pause: hit any key to continue... 13 808 0.2 1.77 1.06 -6.49e-012 1 Pause: hit any key to continue... 14 880 0.2 1.77 1.06 -2.59e-011 1 Pause: hit any key to continue... 15 944 0.2 1.77 1.06 -1.04e-010 1 Pause: hit any key to continue... 16 1003 0.2 1.77 1.06 -4.15e-010 1 Pause: hit any key to continue... 17 1073 0.2 1.77 1.06 -1.66e-009 1 Pause: hit any key to continue... 18 1131 0.2 1.77 1.06 -6.64e-009 1 Pause: hit any key to continue... 19 1193 0.2 1.769 1.06 -2.66e-008 1 Pause: hit any key to continue... 20 1240 0.2 1.769 1.06 -1.06e-007 1 Pause: hit any key to continue... 21 1299 0.200002 1.767 1.06 -4.25e-007 1 Pause: hit any key to continue... 22 1354 0.200007 1.764 1.06 -1.7e-006 1 Pause: hit any key to continue... 23 1408 0.200029 1.759 1.06 -6.82e-006 1 Pause: hit any key to continue... 24 1461 0.200115 1.748 1.05 -2.73e-005 1 Pause: hit any key to continue... 25 1510 0.200462 1.726 1.05 -0.00011 1 Pause: hit any key to continue... 26 1559 0.201858 1.682 1.04 -0.000444 1 Pause: hit any key to continue... 27 1606 0.207495 1.595 1.02 -0.00182 1 Pause: hit any key to continue... 28 1653 0.230524 1.423 0.988 -0.00759 1 Pause: hit any key to continue... 29 1701 0.326949 1.088 0.919 -0.0332 1 Pause: hit any key to continue... 30 1749 0.758103 0.4564 0.785 -0.162 1 Pause: hit any key to continue... 31 1798 1.39346 0 0.223 -1.05 1 Pause: hit any key to continue... 32 1848 1.39346 0 0.35 0 1 Pause: hit any key to continue... xo = 0.8229 0.9114 Ot = 1.3935 nS = 1848 lambda = 1.5945 1.8466 [xo,Ot,nS,lambda]=grg('test14','restr14',[-1 -2],2,1,'gtest14','grestr14',[-5 -5],[5 5]) Pause: hit any key to continue... max Directional Call Iter F-count f(x) constraint Step-size derivative 1 55 13 0 0.05 -324 0 Pause: hit any key to continue... 2 110 3.4 0 0.025 -256 0 Pause: hit any key to continue... 3 165 1.85312 0 0.0141 -64 0 Pause: hit any key to continue... 4 220 1.59673 0 0.00404 -33.1 0 Pause: hit any key to continue... 5 275 1.48984 0 0.00195 -27.9 0 Pause: hit any key to continue... 6 330 1.44046 0 0.000966 -25.8 0 Pause: hit any key to continue... 7 385 1.41668 0 0.000482 -24.8 0 Pause: hit any key to continue... 8 440 1.405 0 0.00024 -24.3 0 Pause: hit any key to continue... 9 495 1.39922 0 0.00012 -24.1 0 Pause: hit any key to continue... 10 550 1.39634 0 6.01e-005 -24 0 Pause: hit any key to continue... 11 604 1.39346 6.904e-007 6.01e-005 -23.9 0 Pause: hit any key to continue... 12 657 1.39346 6.904e-007 0.1 -23.9 1 Pause: hit any key to continue... 13 740 1.39346 6.904e-007 4.35e-017 -3.41 1 Pause: hit any key to continue... 14 823 1.39346 6.904e-007 2.18e-017 -3.41 1 Pause: hit any key to continue... 15 922 1.39346 6.904e-007 3.32e-022 -3.41 1 Pause: hit any key to continue... 16 1005 1.39346 6.904e-007 2.18e-017 -3.41 1 Pause: hit any key to continue... 17 1089 1.39346 6.904e-007 1.09e-017 -3.41 1 Pause: hit any key to continue... 18 1188 1.39346 6.904e-007 3.32e-022 -3.41 1 Pause: hit any key to continue... 19 1270 1.39346 6.904e-007 4.35e-017 -3.41 1 Pause: hit any key to continue... 20 1352 1.39346 6.904e-007 4.35e-017 -3.41 1 Pause: hit any key to continue... 21 1435 1.39346 6.904e-007 2.18e-017 -3.41 1 Pause: hit any key to continue... 22 1518 1.39346 6.904e-007 2.18e-017 -3.41 1 Pause: hit any key to continue... 23 1617 1.39346 6.904e-007 3.32e-022 -3.41 1 Pause: hit any key to continue... 24 1700 1.39346 6.904e-007 2.18e-017 -3.41 1 Pause: hit any key to continue... 25 1783 1.39346 6.904e-007 2.18e-017 -3.41 1 Pause: hit any key to continue... 26 1882 1.39346 6.904e-007 3.32e-022 -3.41 1 Pause: hit any key to continue... 27 1963 1.39346 6.904e-007 8.71e-017 -3.41 1 Pause: hit any key to continue... 28 2045 1.39346 6.904e-007 4.35e-017 -3.41 1 Pause: hit any key to continue... 29 2130 1.39346 6.904e-007 5.44e-018 -3.41 1 Pause: hit any key to continue... 30 2216 1.39346 6.904e-007 2.72e-018 -3.41 1 Pause: hit any key to continue... 31 2303 1.39346 6.904e-007 1.36e-018 -3.41 1 Pause: hit any key to continue... 32 2391 1.39346 6.904e-007 6.8e-019 -3.41 1 Pause: hit any key to continue... 33 2480 1.39346 6.904e-007 3.4e-019 -3.41 1 Pause: hit any key to continue... 34 2571 1.39346 6.904e-007 8.5e-020 -3.41 1 Pause: hit any key to continue... 35 2663 1.39346 6.904e-007 4.25e-020 -3.41 1 Pause: hit any key to continue... 36 2756 1.39346 6.904e-007 2.13e-020 -3.41 1 Pause: hit any key to continue... 37 2850 1.39346 6.904e-007 1.06e-020 -3.41 1 Pause: hit any key to continue... 38 2945 1.39346 6.904e-007 5.31e-021 -3.41 1 Pause: hit any key to continue... 39 3044 1.39346 6.904e-007 3.32e-022 -3.41 1 Pause: hit any key to continue... 40 3127 1.39346 6.904e-007 2.18e-017 -3.41 1 Pause: hit any key to continue... 41 3212 1.39346 6.904e-007 5.44e-018 -3.41 1 Pause: hit any key to continue... 42 3302 1.39346 6.904e-007 1.7e-019 -3.41 1 Pause: hit any key to continue... 43 3395 1.39346 6.904e-007 2.13e-020 -3.41 1 Pause: hit any key to continue... 44 3492 1.39346 6.904e-007 1.33e-021 -3.41 1 Pause: hit any key to continue... 45 3590 1.39346 6.904e-007 6.64e-022 -3.41 1 Pause: hit any key to continue... 46 3689 1.39346 6.904e-007 3.32e-022 -3.41 1 Pause: hit any key to continue... 47 3788 1.39346 6.904e-007 3.32e-022 -3.41 1 Pause: hit any key to continue... 48 3869 1.39346 6.904e-007 8.71e-017 -3.41 1 Pause: hit any key to continue... 49 3968 1.39346 6.904e-007 3.32e-022 -3.41 1 Pause: hit any key to continue... 50 4051 1.39346 6.904e-007 2.18e-017 -3.41 1 Pause: hit any key to continue... Warning GRG: reached maximum number of iterations! xo = 0.8229 0.9114 Ot = 1.3935 nS = 4051 lambda = 1.5945 1.8466 [xo,Ot,nS,lambda]=grg('test14','restr14',[-1 -2],3,1,'gtest14','grestr14',[-5 -5],[5 5]) Pause: hit any key to continue... max Directional Call Iter F-count f(x) constraint Step-size derivative 1 48 793.637 261.2 1.25 -13 0 Pause: hit any key to continue... 2 86 6.07877 9.109 105 -7.57 0 Pause: hit any key to continue... 3 131 3.14087 0 3.99 -1.47 0 Pause: hit any key to continue... 4 183 1.39346 0 0.0581 -25.1 0 Pause: hit any key to continue... 5 233 1.39346 0 0.35 0 0 Pause: hit any key to continue... xo = 0.8229 0.9114 Ot = 1.3935 nS = 233 lambda = 1.5945 1.8466 type grg function [xo,Ot,nS,lambda]=grg(S,Res,x0,ixI,ip,GrS,GrR,Lb,Ub,problem,almax,opt) % Constrained optimization using GRG to solve problems like: % % min S(x) subject to: g(x) <= 0, h(x) = 0 (nonlinear constraints) % x Lb <= x <= Ub % % [xo,Ot,nS,lambda]=grg(S,Res,x0,iXI,ip,GrS,GrR,Lb,Ub,problem,almax,opt) % % S: objective function % x0: initial point % ixI: index vector of independent variables (default: slack variables + first ones) % ip: (0) no plot (default), (>0) plot figure ip with pause, (<0) plot figure ip % Res: contraint function returning [g(x),h(x)] % GrS: gradient vector function of S(x) % GrR: gradient matrix function of Res(x), returning [Dg(x),Dh(x)] % Lb, Ub: lower and upper bound vectors [used also to plot (plot default = x0*(1+/-2))] % problem: (-1): minimum (default), (1): maximum % almax: maximum stepsize (default = 10) % opt: options vector set using optimset('grg') % xo: optimal point % Ot: optimal value of S % nS: number of objective function evaluations % lambda: Lagrange multipliers % Copyright (c) 2005 by LASIM-DEQUI-UFRGS % $Revision: 1.0 $ $Date: 2005/09/09 08:00:00 $ % Euclides Almeida Neto e Argimiro R. Secchi (arge@enq.ufrgs.br) % [xo,Ot,nS,lambda]=grg('test1','restr1',[-1.2 1],[1 2],-1,'gtest1','grestr1',[-5 -6],[5 6]) % [xo,Ot,nS,lambda]=grg('test14','restr14',[-1 -2],1,-1,'gtest14','grestr14',[-6 -6],[6 6]) % [xo,Ot,nS,lambda]=grg('test14','restr14eq',[-1 -2 0],1,0,'gtest14eq','grestr14eq',[-6 -6 0],[6 6 inf]) % [xo,Ot,nS,lambda]=grg('test15','restr15',[-1 -2],1,-1,'gtest15','grestr15',[-6 -6],[6 6]) % [xo,Ot,nS,lambda]=grg('test16','restr16',[0.05 0.5],[1 2],-1,'gtest16','grestr16',[0 0],[1 2]) % [xo,Ot,nS,lambda]=grg('test17','restr17',[-1 1],[1 3],-1,'gtest17','grestr17',[-5 -6],[5 6]) % [xo,Ot,nS,lambda]=grg('test23','restr23',[2 2 2],1,0,'gtest23','grestr23',[],[],[],1) % [xo,Ot,nS,lambda]=grg('test24','restr24',[0 0],[1 2],-1,'gtest24','grestr24',[-5 -5],[5 5]) % [xo,Ot,nS,lambda]=grg('test24','restr24eq',[0 0 0 0],[1 2],0,'gtest24eq','grestr24eq',[-5 -5 0 0],[5 5 Inf Inf]) % If just 'defaults' passed in, return the default options in xo if nargin==1 & nargout <= 1 & isequal(S,'defaults') defaultopt = optimset('display','final','TolX',1e-6,'TolFun',1e-6,'TolCon',1e-6,... 'DerivativeCheck','off','Diagnostics','off','GradObj','on','GradConstr','on',... 'MaxIter','50*(1+4*~(ip>0))'); xo = defaultopt; return end global X_PLOT h1 h2 % Definitions INIT = 1; OUTERLOOP = 2; INNERLOOP = 3; NEWTONLOOP = 4; END_INNERLOOP = 5; FINAL_SOLUTION = 6; TRUE = 1; FALSE = 0; % Verifying Arguments if nargin < 3, error('GRG requires three input arguments!'); end x0=x0(:); if nargin < 4 | isempty(ixI), ixI=[]; else ixI=ixI(:); end if nargin < 5 | isempty(ip), ip=0; end if nargin < 6 | isempty(GrS), gradS=0; else gradS=1; end if nargin < 7 | isempty(GrR), gradR=0; else gradR=1; end if nargin < 8 | isempty(Lb), Lbx=-Inf*ones(length(x0),1); Lb=-x0-~x0; else Lbx=Lb(:); end if nargin < 9 | isempty(Ub), Ubx=Inf*ones(length(x0),1); Ub=2*x0+~x0; else Ubx=Ub(:); end if nargin < 10 | isempty(problem), problem=-1; end if nargin < 11 | isempty(almax), almax = 1e10; end if nargin < 12 | isempty(opt), opt = optimset('grg'); end if isfield(opt,'Current_nS'), nS = opt.Current_nS; rec = opt.Recur; Par.it = opt.Iterations; else nS = 0; % Number of Objective Function Evaluations, first call rec = 0; % Number of recurrent calls Par.it = 0; % Iteration of outer loop h1 = []; h2 = []; end % -------------------------------------- % Initialization of the GRG Algoritm % -------------------------------------- mxit = eval(optimget(opt,'MaxIter')); TolX = optimget(opt,'TolX'); TolFun = optimget(opt,'TolFun'); TolCon = optimget(opt,'TolCon'); alfa = almax; Par.tol = TolFun; Par.gd = -Inf; Par.n2gd = Inf; Par.itn = 0; Par.iti = 0; nx = length(x0); % Number of Variables switch optimget(opt,'Display') case {'iter'} verbosity = 3; case {'final'} verbosity = 2; case {'off'} verbosity = 1; otherwise verbosity = 0; end formatstr = '%5.0f %5.0f %12.6g %12.4g %12.3g %12.3g %d'; if ip & nx == 2 & nS == 0, figure(abs(ip)); [X1,X2]=meshgrid(Lb(1):(Ub(1)-Lb(1))/20:Ub(1),Lb(2):(Ub(2)-Lb(2))/20:Ub(2)); [n1,n2]=size(X1); f=zeros(n1,n2); [ine,eq]=feval(Res,x0); ni=size(ine(:),1); ne=size(eq(:),1); if ni, ri=zeros(n1,n2,ni); end if ne, re=zeros(n1,n2,ne); end for i=1:n1, for j=1:n2, f(i,j)=feval(S,[X1(i,j);X2(i,j)]); [ine,eq]=feval(Res,[X1(i,j);X2(i,j)]); if ni, ri(i,j,:)=ine; end if ne, re(i,j,:)=eq; end end end mxf=max(max(f)); mnf=min(min(f)); df=mnf+(mxf-mnf)*(2.^(([0:10]/10).^2)-1); [v,h]=contour(X1,X2,f,df); hold on; if ni, for i=1:ni, contour(X1,X2,ri(:,:,i),[0 0],'r:'); end end if ne, for i=1:ne, contour(X1,X2,re(:,:,i),[0 0],'r:'); end end % clabel(v,h); % h1=plot(x0(1),x0(2),'ro'); h1=xplot(x0,0); legend(h1,'start point'); if ip > 0, disp('Pause: hit any key to continue...'); pause; end X_PLOT=1; end if isempty(rec), rec = 0; h2 = 0; end if problem == 1, pr='-feval('''; else pr='feval('''; end; fun=[pr S ''',x)']; % ======================================================================================= % 1) Choosing a Feasible point x0 and Selecting n - m Independents Variables, k = 0 % ======================================================================================= Snew = feval(S,x0); Sx = Inf; nS = nS + 1; [g,he]=feval(Res,x0); g=g(:); he=he(:); s0 = -g; % Initial value of "s" if ~isempty(g), h = [he; g+s0]; else h = he; end % ----------------------------------------------- % Dimension of the Optimization Problem. % ----------------------------------------------- nh = size(he,1); % Number of Equality Constraints ng = size(g,1); % Number of Inequality Constraints n = nx + ng; % Number of Variables including Slacks nD = nh + ng; % Number of Dependent Variables nI = n - nD; % Number of Independent Variables dgs = eye(ng); if ~nh | ~ng zer=[]; else zer=zeros(nh,ng); end LB = [Lbx(:); zeros(ng,1)]; UB = [Ubx(:); Inf*ones(ng,1)]; OPT = optimset('fminbnd'); % linesearch OPT = optimset(OPT,'Display','off','TolX',TolX); x = [x0; s0]; % ---------------------------------------- % Verifying Feasibility of the point x0. % ---------------------------------------- nf = find(abs(h)>Par.tol); if isempty(nf) & LB-x < Par.tol & x-UB < Par.tol Par.feasible = TRUE; else Par.feasible = FALSE; end % ---------------------------------------------------- % Looking for Indexes of the Independent Variables % ---------------------------------------------------- iv=length(ixI); if iv < nI, if iv > 0, error('Invalid number of independent variables!'); end if ng, idx=[nx+1:nx+ng 1:nx]'; else idx=[1:nx]'; end ixI=idx(1:nI); end ixD = []; for i = 1: n xx = find(ixI == i); zz = isempty(xx); if (zz == 1) ixD = [ixD; i]; end end % ---------------------------------------------------- % Analyzing the Feasibility of the Initial Point - x0 % ---------------------------------------------------- if verbosity > 2, Par.alfa = alfa; Par.status = INIT; ShowData(x0,Snew,he,g,nS,Par,rec); end if verbosity > 1 & nS == 1, header = sprintf(['\n max Directional Call\n',... ' Iter F-count f(x) constraint Step-size derivative']); disp(header); end xi = x0; s = s0; itx = 0; ineq=[s; x-LB; UB-x]; nsv=find(ineq<0); mg = sum(max(abs(he))) + sum(max(abs(ineq(nsv)))); % sum: to zero empty vector % =-=-=-=-=-=-=-=-=-=-=-=-=- Outer Loop =-=-=-=-=-=-=-=-=-=-=-=-=- while (Par.n2gd > 1 & Par.it < mxit) % Outer Loop Par.it = Par.it + 1; Sx = Snew; if verbosity > 2, Par.status = OUTERLOOP; ShowData(xi,Snew,he,g,nS,Par,rec); end if ip & nx == 2 & Par.it > 1, xplot(xi,ip); end % ================================================================================ % 2) Linearization of the Objective Function and Constraints around the point xk % ================================================================================ % --------------------------- % Evaluation of dS(x) % --------------------------- dS = feval(GrS,xi); dS = dS(:)'; % ------------------------------ % Evaluation of dh(x) and dg(x) % ------------------------------ [dg,dhe] = feval(GrR,xi); if ~isempty(g), h = [he; g+s]; else h = he; end nf = find(abs(h)>Par.tol); if isempty(nf) & LB-x < Par.tol & x-UB < Par.tol Par.feasible = TRUE; else; Par.feasible = FALSE; end dh = [dhe zer; dg dgs]; b = dh * x - h; % -------------------------------------------------------------------------- % Calculation of dDS(x), dIS(x), B and C (D - Dependents, I - Independents) % -------------------------------------------------------------------------- dS1 = [dS zeros(1,ng)]; dDS = dS1(:,ixD)'; dDh = dh(:,ixD); xD = x(ixD); dIS = dS1(:,ixI)'; dIh = dh(:,ixI); xI = x(ixI); B = dDh; C = dIh; % ================================================================================ % 3) Calculation of Reduced Gradient - gR % ================================================================================ iB = inv(B); lambda = -iB'*dDS; gR = dIS + C'*lambda; % ================================================================================ % 4) Calculation of Search Direction in the Independent Variables Space % ================================================================================ dI = -gR; % dI = []; % for i = 1: nI % if ((abs(xI(i) - UB(ixI(i))) < Par.tol & gR(i) < 0) | (abs(xI(i) - LB(ixI(i))) < Par.tol & gR(i) > 0)) % colocar Inclusive vars de folga % dI = [dI; 0]; % else % dI = [dI; -gR(i)]; % end % end % ================================================================================ % 5) Calculation of Search Direction in the Dependent Variables Space % ================================================================================ dD = -iB*C*dI; % ================================================================================ % 6) Linesearch % ================================================================================ dDI = [xI dI ixI]; dDD = [xD dD ixD]; dDI = [dDD; dDI]; dDI = sortrows(dDI,3); x = dDI(:,1); xi = x(1:nx); d = dDI(:,2); di = d(1:nx); [alfa,Sal,fl,out] = fminbnd('ObjFunSa',1e-10,almax,OPT,S,xi,di); nS = nS + out.iterations; for i = 1:n if d(i) ~= 0 alfa_LB = (LB(i) - x(i))/d(i); alfa_UB = (UB(i) - x(i))/d(i); alfab = max([alfa_LB, alfa_UB]); if (alfab <= Par.tol) % error('change independent variables'); k=find(ixI==i); if isempty(k) if alfab, ixI = [i; ixI]; rec = rec + 1; if verbosity > 1, CurrOutput = sprintf(formatstr,Par.it,nS,Snew,mg,alfa,Par.gd,rec); disp(CurrOutput); end if rec > n disp('maximum changes of independent variables'); xo = xi; Ot = Snew; return; end if ip & nx == 2, h2=plot(x(1),x(2),'kh'); if rec==1, legend([h1,h2],'start point','new ind. var.'); end end opt.Current_nS = nS; opt.Recur = rec; opt.Iterations = Par.it; [xo,Ot,nS,lambda]=grg(S,Res,xi,ixI(1:nI),ip,GrS,GrR,Lbx,Ubx,problem,almax,opt); return; end else if alfab, dI(k) = -d(i); gR(k) = -gR(k); % change signal to avoid early exit alfa = min([alfa, -alfab]); else dI(k) = 0; end end else alfa = min([alfa, alfab]); end end end % ================================================================================ % 7) Calculation of k % ================================================================================ % +=+=+=+=+=+=+=+=+=+=+=+=+= Inner Loop +=+=+=+=+=+=+=+=+=+=+=+=+= Par.gd = gR'*dI; % ================================================================================================== % 9) If xk+1 is Feasible and S(xk+1) <= S(xk), then k <- k + 1, and If the convergence criterion % is not satisfied (go to step 2) else END. % ================================================================================================== Par.iti = 0; feas = Par.feasible; xo = xi; while ~Par.iti | ((feas & (~Par.feasible | Snew-Sx > 0*Par.tol)) & Par.iti < mxit) % Inner Loop Par.iti = Par.iti + 1; if verbosity > 2, Par.status = INNERLOOP; Par.alfa = alfa; ShowData(x,Snew,he,g,nS,Par,rec); end xD1 = xD + alfa * dD; xI1 = xI + alfa * dI; % ====================================================================================================== % 8) Using Newton method to push Dependent Variables to the Feasible Region h(xk+1) = 0 (correction), % utilizing xDk+1 as Initial Estimation. % ====================================================================================================== alfa = alfa/2; % ------------------------------------------- % Newton Method to the Algebraic Equations % ------------------------------------------- sx = [xI1 ixI; xD1 ixD]; xID = sortrows(sx,2); x = xID(:,1); xi = x(1:nx); if n > nx s = x(nx+1:n); end error(1:3) = inf; Par.itn = 0; while (error(1) > Par.tol & Par.itn < mxit) % Newton Step Loop Par.itn = Par.itn + 1; [g,he] = feval(Res,xi); g=g(:); he=he(:); [dg,dhe] = feval(GrR,xi); if ~isempty(g), h = [he; g+s]; else h = he; end dh = [dhe zer; dg dgs]; dDh = dh(:,ixD); dx = -inv(dDh)*h; xD1 = xD1 + dx; error(3) = error(2); error(2) = error(1); error(1) = norm(dx)/(norm(xD1) + eps) + norm(h); if error(3) < error(2) & error(2) < error(1), % to prevent newton divergence break; end sx = [xI1 ixI; xD1 ixD]; xID = sortrows(sx,2); x = xID(:,1); xi = x(1:nx); if n > nx s = x(nx+1:n); end end % Newton Step Loop % -+-+-+-+-+-+-+-+-+-+-+-+-+ Newton Step -+-+-+-+-+-+-+-+-+-+-+-+- nf = find(abs(h)>Par.tol); if isempty(nf) & LB-x < Par.tol & x-UB < Par.tol Par.feasible = TRUE; else Par.feasible = FALSE; end Snew = feval(S,xi); nS = nS + 1; [g,he]=feval(Res,xi); g=g(:); he=he(:); % ================================================================================================= % 10) If xk+1 is Feasible and S(xk+1) > S(xk) or Newton Method din't converged, % then reduce alfak and (go to step 7). % ================================================================================================= if verbosity > 2, Par.status = NEWTONLOOP; ShowData(xi,Snew,he,g,nS,Par,rec); end end % Inner Loop % +=+=+=+=+=+=+=+=+=+=+=+=+= Inner Loop +=+=+=+=+=+=+=+=+=+=+=+=+= if ~Par.feasible & ng, x=[xi;-g]; end ineq=[s; x-LB; UB-x]; nsv=find(ineq<0); mg = sum(max(abs(he))) + sum(max(abs(ineq(nsv)))); % sum: to zero empty vector Par.n2gd = abs(Par.gd); if verbosity > 1, CurrOutput = sprintf(formatstr,Par.it,nS,Snew,mg,alfa,Par.gd,rec); disp(CurrOutput); if verbosity > 2, Par.status = END_INNERLOOP; ShowData(xi,Snew,he,g,nS,Par,rec); end end Par.n2gd = Par.n2gd/TolFun + norm(xo-xi)/TolX + mg/TolCon; % convergence criteria end % Outer Loop % =-=-=-=-=-=-=-=-=-=-=-=-=- Outer Loop =-=-=-=-=-=-=-=-=-=-=-=-=- % Showing Optimun Solution xo = xi; Ot = Snew; if verbosity > 2, Par.status = FINAL_SOLUTION; ShowData(xi,Snew,he,g,nS,Par,rec); end if ip & nx == 2, xplot(xi,ip); X_PLOT=0; end if Par.it < mxit, if ip & nx == 2, figure(abs(ip)); h3=plot(xo(1),xo(2),'r*'); if h2, legend([h1,h2,h3],'start point','new ind. var.','optimum'); else legend([h1,h3],'start point','optimum'); end end else disp('Warning GRG: reached maximum number of iterations!'); end edit grg disp(sprintf('%d',feature('SessionTool'))) 0 ; dbstatus dbstack dbstack ; disp(which('grg')); V:\Cursos\Pos\Otimiza\Aulas\grg.m mdbstatus 'v:\cursos\pos\otimiza\aulas\grg.m' ;