Home > source > @penlab > eqconstr_min.m

# eqconstr_min

## PURPOSE

equality constrained minimization

## SYNOPSIS

function [nFlag,rResults]=eqconstr_min(obj)

## DESCRIPTION

     equality constrained minimization

     obj() etc. should work with parameters x and UEQ (equality lagrangian
     multipliers). If UEQ!=0, the function returns lagrangian
     (e.i., augmented lagrangian for inequalities + lagrangian for equalities)
     as the first return value and function value/gradient of equality constraints
     as the second one; obj_hes() has just one return argument (hessian)
     expects an open nl-file?, N, N_EQUAL, ... global variables

     the equality constrainted problem
     update different LS procedures

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- penlab PenLab (Pennon Laboratory for Matlab, previously PennonM)

## SOURCE CODE

    0001 function [nFlag,rResults]=eqconstr_min(obj)
    0002 % equality constrained minimization
    0003 %
    0004 % obj() etc. should work with parameters x and UEQ (equality lagrangian
    0005 % multipliers). If UEQ!=0, the function returns lagrangian
    0006 % (e.i., augmented lagrangian for inequalities + lagrangian for equalities)
    0007 % as the first return value and function value/gradient of equality constraints
    0008 % as the second one; obj_hes() has just one return argument (hessian)
    0009 % expects an open nl-file?, N, N_EQUAL, ... global variables
    0010 %
    0011 % the equality constrainted problem
    0012 % update different LS procedures
    0013
    0014
    0015 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0016
    0017   starttime=cputime;
    0018   % reset default parameters if necessary
    0019   MAX_MITER = obj.allopts.max_inner_iter;
    0020   ALPHA = obj.allopts.inner_stop_limit;
    0021   TOL_DIR = obj.allopts.eq_dir_stop_limit;
    0022   solver = obj.allopts.eq_solver;
    0023   linesearch = obj.allopts.eq_linesearch;
    0024   solver_warn_max = obj.allopts.eq_solver_warn_max;
    0025   ls_short_max = obj.allopts.ls_short_max;
    0026   recover_strategy = obj.allopts.min_recover_strategy;
    0027   recover_max = obj.allopts.min_recover_max;
    0028
    0029   miter=0;
    0030   solver_warn=0;
    0031   ls_short=0;
    0032   recover_no=0;
    0033   tol_dir_orig=TOL_DIR;
    0034
    0035   %x=x0;       starting point is the current point
    0036   %ueq = ueq0;
    0037   %[fx,eqx]=fnc.obj(x,ueq);
    0038   %[gradx, gradeqx]=fnc.obj_grad(x,ueq);
    0039   obj.eval_alx();
    0040   obj.eval_aldx();
    0041   fx_start = obj.ALx;
    0042   rNormG = norm(obj.ALdx);
    0043   rNormG_start = rNormG;
    0044   infeas_eq = norm(obj.eqx);
    0045   infeas_eq_start = infeas_eq;
    0046
    0047   obj.print(3,Inf,'object(x_%3i) = %24.16E',miter,obj.ALx);
    0048   obj.print(3,Inf,'||grad(x)||_2 = %24.16E',rNormG);
    0049   obj.print(3,Inf,'infeas_eq     = %24.16E',infeas_eq);
    0050   obj.print(3,4,'     ----');
    0051   obj.print(4,Inf,'        --- start of inner iter ---');
    0052   obj.print(4,Inf,' ');
    0053
    0054   nFlag=1; % max it limit reached
    0055   rResults=[];
    0056   obj.initer_last=0;
    0057   obj.lsiter_last=0;
    0058   obj.stats_time_fact_last=0;
    0059
    0060   % init LS
    0061   if (linesearch==9)
    0062     error('not adjusted for PenLab yet');
    0063     %lseq_inoc(fnc);
    0064     obj.lseq_inoc();
    0065   end
    0066
    0067     while (miter < MAX_MITER)
    0068
    0069       %hessx=fnc.obj_hess(x,ueq);
    0070       obj.eval_alddx();
    0071       %hessx=obj.ALddx;
    0072
    0073       %%%%%%%% Newton solver %%%%%%%%%
    0074       switch solver
    0075       case 0  % LDL factorization
    0076         [x_dir,ueq_dir,nFlagSol] = obj.solvekkt_ldl(obj.ALddx,obj.eqdx,-obj.ALdx,-obj.eqx);
    0077       case 1  % LU factorization
    0078         [x_dir,ueq_dir,nFlagSol] = obj.solvekkt_lu(obj.ALddx,obj.eqdx,-obj.ALdx,-obj.eqx);
    0079       case 2  % MA57 factorization
    0080         [x_dir,ueq_dir,nFlagSol] = obj.solvekkt_ma57(obj.ALddx,obj.eqdx,-obj.ALdx,-obj.eqx);
    0081       case 3  % Luksan 3
    0082         [x_dir,ueq_dir,nFlagSol] = obj.solvekkt_luksan3(obj.ALddx,obj.eqdx,-obj.ALdx,-obj.eqx,TOL_DIR,20+obj.Nx+obj.NYnnz+obj.Neq);
    0083       case 4  % Luksan 4
    0084         [x_dir,ueq_dir,nFlagSol] = obj.solvekkt_luksan4(obj.ALddx,obj.eqdx,-obj.ALdx,-obj.eqx,TOL_DIR,20+2*(obj.Nx+obj.NYnnz-obj.Neq));
    0085       case 5  % Schur complement + CGM
    0086         [x_dir,ueq_dir,nFlagSol] = obj.solvekkt_schurcgm(obj.ALddx,obj.eqdx,-obj.ALdx,-obj.eqx,TOL_DIR,20+2*obj.Neq);
    0087       case 6  % Projected CGM onto nullspace (implicit)
    0088         [x_dir,ueq_dir,nFlagSol] = obj.solvekkt_projcgm(obj.ALddx,obj.eqdx,-obj.ALdx,-obj.eqx,TOL_DIR,20+2*(obj.Nx+obj.NYnnz-obj.Neq));
    0089       case 7  % Projected CGM onto nullspace (implicit, LDL)
    0090         [x_dir,ueq_dir,nFlagSol] = obj.solvekkt_projcgm2(obj.ALddx,obj.eqdx,-obj.ALdx,-obj.eqx,TOL_DIR,20+2*(obj.Nx+obj.NYnnz-obj.Neq));
    0091       case 8  % CGM on the nullspace (implicit Z based on the structure of FMO)
    0092         [x_dir,ueq_dir,nFlagSol] = obj.solvekkt_nscgm(obj.ALddx,obj.eqdx,-obj.ALdx,-obj.eqx,TOL_DIR,20+2*(obj.Nx+obj.NYnnz-obj.Neq));
    0093       case 103  % Luksan 3 with a special stop crit
    0094         stop_test=@obj.lseq_inoc;
    0095         [x_dir,ueq_dir,nFlagSol] = obj.solvekkt_luksan3e(obj.ALddx,obj.eqdx,-obj.ALdx,-obj.eqx,stop_test,20+obj.Nx+obj.NYnnz+obj.Neq);
    0096       case 104  % Luksan 4 & LS stop crit
    0097         stop_test=@obj.lseq_inoc;
    0098         [x_dir,ueq_dir,nFlagSol] = obj.solvekkt_luksan4e(obj.ALddx,obj.eqdx,-obj.ALdx,-obj.eqx,stop_test,20+2*(obj.Nx+obj.NYnnz-obj.Neq));
    0099       %case   %
    0100       otherwise
    0101         obj.print(1,Inf,'ERR @ eqconstr_min(): Sorry, no such solver, terminating...');
    0102         nFlag=100;
    0103         break;
    0104       end
    0105
    0106       if (nFlagSol>0)  % solver fatal error
    0107         nFlag=2; % solver failed, cannot continue
    0108         obj.print(3,Inf,'FAILURE: KKT solver, cannot continue (%i)',nFlagSol);
    0109         break;
    0110       elseif (nFlagSol<0)  % solver error, however, hopefully recoverable
    0111         solver_warn = solver_warn+1;
    0112         if (solver_warn>solver_warn_max)
    0113           nFlag=2;
    0114           obj.print(3,Inf,'FAILURE: KKT solver, cannot continue (flag %i or similar happened %i times)',nFlagSol,solver_warn);
    0115           break;
    0116         else
    0117           obj.print(4,Inf,'Step accepted as a trial one (remaining: %d).\n',solver_warn_max-solver_warn);
    0118         end
    0119       else % solver OK
    0120         solver_warn=0;
    0121       end
    0122
    0123       %%%%%%%% linesearch %%%%%%%%%
    0124       % get the maximal step length
    0125       %rMaxStep=obj.get_maxstep(x_dir);   % TODO ... no barriers so no problem
    0126       rMaxStep=1;
    0127       nMaxSteps=ceil(-log(rMaxStep)/log(2));
    0128       obj.print(5,Inf,'LSEQ: max step %.1e, %i steps',rMaxStep,nMaxSteps);
    0129
    0130       switch linesearch
    0131       case 0  % Do nothing, leave original data
    0132      % such as x, grad_x (useful for TR, ...)
    0133       case 1  % Do full steps, no linesearch at all
    0134         [rRelStep, nFlagLS] = obj.lseq_fullstep(x_dir,ueq_dir);
    0135       %case 2  % Armijo linesearch
    0136       %  [rRelStep, nFlagLS] = ls_armijo(dir);
    0137       case 3  % Pennlp/Pennon equality linesearch
    0138         [rRelStep, nFlagLS] = obj.lseq_pen(x_dir,ueq_dir,miter);
    0139       case 4  % Filter linesearch, first version, Biegler/Wachter
    0140         [rRelStep, nFlagLS] = obj.lseq_filter(x_dir,ueq_dir,miter);
    0141       case 5  % Nocedal Knitro LS, simple, first version
    0142         [rRelStep, nFlagLS] = obj.lseq_noc(x_dir,ueq_dir,obj.ALddx*x_dir,miter);
    0143       case 6  % Nocedal Flexible Merit LS
    0144         [rRelStep, nFlagLS] = obj.lseq_flex(x_dir,ueq_dir,obj.ALddx*x_dir,miter);
    0145       case 7  % Nocedal Knitro LS, simple, 2nd version
    0146         [rRelStep, nFlagLS] = obj.lseq_noc2(x_dir,ueq_dir,obj.ALddx*x_dir,miter);
    0147       case 8  % Nocedal Knitro LS, simple, 3nd version
    0148         [rRelStep, nFlagLS] = obj.lseq_noc3(x_dir,ueq_dir,obj.ALddx*x_dir,miter);
    0149       case 9  % Nocedal Knitro inexact SQP LS
    0150         [nFlagLS, rRelStep] = obj.lseq_inoc(x_dir,ueq_dir,obj.ALddx*x_dir);
    0151       %case   df%
    0152       otherwise
    0153         obj.print(1,Inf,'ERR @ eqconstr_min(): Sorry, no such LS, terminating...');
    0154         nFlag=100;
    0155         break;
    0156       end
    0157
    0158       if (nFlagLS>0 && solver_warn>0)
    0159         nFlag=2; % LS failed because of solver, cannot continue
    0160         obj.print(3,Inf,'FAILURE: Linesearch (because of the KKT solver), cannot continue');
    0161         %break;
    0162       elseif (nFlagLS>0)
    0163         nFlag=3; % LS failed, cannot continue
    0164         obj.print(3,Inf,'FAILURE: Linesearch, cannot continue');
    0165         %break;
    0166       elseif (nFlagLS<0)
    0167         ls_short = ls_short+1;
    0168         if (ls_short > ls_short_max)
    0169           nFlag=3; % LS failed, cannot continue
    0170           obj.print(3,Inf,'FAILURE: Linesearch (flag %d, wawrnings %d times), cannot continue',nFlagLS,ls_short);
    0171           %break;
    0172         else
    0173           obj.print(4,Inf,'Short LS step accepted (remaining %d)',ls_short_max-ls_short);
    0174         end
    0175       else  % LS step OK
    0176         ls_short=0;
    0177       end
    0178
    0179       % return original settings
    0180       if (recover_no>0)
    0181         %TOL_DIR = check_opt(PEN_OPT, 'eq_dir_stop_limit', TOL_DIR_def);
    0182         TOL_DIR=tol_dir_orig/10;
    0183         solver = obj.allopts.eq_solver;
    0184         if (nFlag<=0)  % strategy was successful
    0185           recover_no=0;
    0186         end
    0187       end
    0188
    0189       % recover from LS/solver failure?
    0190       if (nFlag>1)  % not for iterarion limit
    0191         switch recover_strategy
    0192         case 0  % do nothing - fail
    0193           break;
    0194
    0195         case 1  % increase precision of the iterative solver or switch to LDL
    0196           if (recover_no>=recover_max)  % failed to recover
    0197             obj.print(3,Inf,'FAILURE: recovery strategy - no attemps left %i out of %i',recover_no, recover_max);
    0198             break;
    0199           end
    0200           if (solver<3)  % applicable only for iterative methods
    0201             break;
    0202           end
    0203           if (solver_warn>0 || recover_no>=recover_max-1)   % solver in troubles -> do not increase precision -> switch solver
    0204             obj.print(3,Inf,'\nRecover strategy - switching from iterative methods to LDL');
    0205             solver=0;
    0206           else
    0207             TOL_DIR=TOL_DIR/(100^(recover_no+1));
    0208           end
    0209
    0210         case 2  % switch to LDL
    0211           if (recover_no>0)
    0212             obj.print(3,Inf,'FAILURE: recover strategy to switch to LDL didn''t help');
    0213             break;
    0214           end
    0215           if (solver==0)
    0216             break;
    0217           end
    0218           solver=0;
    0219
    0220         otherwise
    0221           obj.print(3,Inf,'WARNING @ eqconstr_min(): Sorry, no such min_recover_strategy');
    0222           break;
    0223         end
    0224         % try to recover with the new settings
    0225         recover_no = recover_no+1;
    0226         obj.print(3,Inf,'\nTry to recover from failure (nFlag=%i) using strategy %i',nFlag,recover_strategy);
    0227         nFlag=0;
    0228         continue;
    0229       end
    0230
    0231
    0232       %%%%%%%% %%%%%%%%%
    0233       rNormG = norm(obj.ALdx);
    0234       infeas_eq = norm(obj.eqx);
    0235       miter=miter+1;
    0236
    0237       obj.print(4,Inf,' ');
    0238       obj.print(3,Inf,'object(x_%3i) = %24.16E',miter,obj.ALx);
    0239       obj.print(3,Inf,'||grad(x)||_2 = %24.16E',rNormG);
    0240       obj.print(3,Inf,'infeas_eq     = %24.16E',infeas_eq);
    0241       obj.print(3,4,'     ----');
    0242       obj.print(4,Inf,'        --- end of %3i in. iter ---\n',miter);
    0243
    0244       %%%%%%%% stopping criterion %%%%%%%%%
    0245       % just make it simple at the begining
    0246       if (rNormG < ALPHA && infeas_eq < ALPHA)
    0247         nFlag=0;
    0248         break;
    0249       end
    0250
    0251     end % of while
    0252
    0253
    0254   rResults=[obj.ALx,rNormG, rNormG_start];
    0255
    0256   % update stats
    0257   obj.initer = obj.initer+obj.initer_last;
    0258   obj.lsiter = obj.lsiter+obj.lsiter_last;
    0259   obj.miter=obj.miter+miter;
    0260   obj.miter_last=miter;
    0261   obj.stats_time_miter_last = cputime - starttime;
    0262   obj.stats_time_miters = obj.stats_time_miters + obj.stats_time_miter_last;
    0263   obj.stats_time_fact = obj.stats_time_fact + obj.stats_time_fact_last;
    0264
    0265
    0266   return;
    0267