matRad_OptimizerFmincon

0001 classdef matRad_OptimizerFmincon < matRad_Optimizer
0002 % matRad_OptimizerFmincon implements the interface for the fmincon optimizer
0003 % of the MATLAB Optiization toolbox
0004 %
0005 % References
0006 %   -
0007 %
0008 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0009 %
0010 % Copyright 2019 the matRad development team.
0011 %
0012 % This file is part of the matRad project. It is subject to the license
0013 % terms in the LICENSE file found in the top-level directory of this
0014 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0015 % of the matRad project, including this file, may be copied, modified,
0016 % propagated, or distributed except according to the terms contained in the
0017 % LICENSE file.
0018 %
0019 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0020 
0021     
0022     properties
0023         options     %the optimoptions for fmincon
0024         wResult     %last optimization result
0025         resultInfo  %info struct about last results
0026     end
0027     
0028     methods
0029         function obj = matRad_OptimizerFmincon
0030             %matRad_OptimizerFmincon
0031             %   Construct an instance of the fmincon optimizer from the Optimization Toolbox
0032             matRad_cfg = MatRad_Config.instance();
0033             
0034             obj.wResult = [];
0035             obj.resultInfo = [];
0036             
0037             %createDefaultOptimizerOptions Constructs a set of default
0038             %options for the optimizer to use
0039             obj.options = optimoptions('fmincon',...
0040                 'Algorithm','interior-point',...
0041                 'Display','iter-detailed',...
0042                 'SpecifyObjectiveGradient',true,...
0043                 'SpecifyConstraintGradient',true,...
0044                 'AlwaysHonorConstraints', 'bounds',...
0045                 'MaxIterations',matRad_cfg.propOpt.defaultMaxIter,...
0046                 'MaxFunctionEvaluations',3000,...
0047                 'CheckGradients',false,...
0048                 'HessianApproximation',{'lbfgs',6},...
0049                 'UseParallel',true,...
0050                 'Diagnostics','on',...
0051                 'ScaleProblem',true,...
0052                 'PlotFcn',{@optimplotfval,@optimplotx,@optimplotfunccount,@optimplotconstrviolation,@optimplotstepsize,@optimplotfirstorderopt});                    
0053         end
0054                 
0055         function obj = optimize(obj,w0,optiProb,dij,cst)
0056             %optimize Carries Out the optimization
0057             
0058             % obtain lower and upper variable bounds
0059             lb = optiProb.lowerBounds(w0);
0060             ub = optiProb.upperBounds(w0);
0061                         
0062             % Informing user to press q to terminate optimization
0063             %fprintf('\nOptimzation initiating...\n');
0064             %fprintf('Press q to terminate the optimization...\n');
0065             
0066             % Run fmincon.
0067             [obj.wResult,fVal,exitflag,info] = fmincon(@(x) obj.fmincon_objAndGradWrapper(x,optiProb,dij,cst),...
0068                 w0,... % Starting Point
0069                 [],[],... % Linear Constraints we do not explicitly use
0070                 [],[],... % Also no linear inequality constraints
0071                 lb,ub,... % Lower and upper bounds for optimization variable
0072                 @(x) obj.fmincon_nonlconWrapper(x,optiProb,dij,cst),...
0073                 obj.options); % Non linear constraint structure);
0074             
0075             obj.resultInfo = info;
0076             obj.resultInfo.fVal = fVal;
0077             obj.resultInfo.exitflag = exitflag;
0078         end
0079         
0080         function [f, fGrad] = fmincon_objAndGradWrapper(obj,x,optiProb,dij,cst)
0081             f = optiProb.matRad_objectiveFunction(x,dij,cst);
0082             fGrad = optiProb.matRad_objectiveGradient(x,dij,cst);
0083         end
0084         
0085         function [c,cEq,cJacob,cEqJacob] = fmincon_nonlconWrapper(obj,x,optiProb,dij,cst)
0086             %Get the bounds of the constraint
0087             [cl,cu] = optiProb.matRad_getConstraintBounds(cst);
0088                     
0089             %Get finite bounds
0090             clFinIx = isfinite(cl);
0091             cuFinIx = isfinite(cu);
0092             
0093             % Some checks
0094             assert(isequal(size(cl),size(cu)));
0095             assert(all(cl <= cu));
0096             
0097             %For fmincon we need to separate into equalty and inequality
0098             %constraints
0099             isEqConstr = (cl == cu);
0100             eqIx = isEqConstr;
0101             ineqIx = ~isEqConstr;
0102             
0103             %Obtain all constraint functions and derivatives
0104             cVals = optiProb.matRad_constraintFunctions(x,dij,cst);
0105             cJacob = optiProb.matRad_constraintJacobian(x,dij,cst);
0106             
0107             %Subselection of equality constraints
0108             cEq = cVals(eqIx & clFinIx); %We can only rely on cl indices here due to the equality index
0109             cEqJacob = cJacob(eqIx & clFinIx,:)';
0110             
0111             %Prepare inequality constraints:
0112             %We need to separate upper and lower bound constraints for
0113             %fmincon
0114             cL = cl(ineqIx & clFinIx) - cVals(ineqIx & clFinIx);
0115             cU = cVals(ineqIx & cuFinIx) - cu(ineqIx & cuFinIx);
0116             cJacobL = -cJacob(ineqIx & clFinIx,:);
0117             cJacobU = cJacob(ineqIx & cuFinIx,:);
0118             
0119             %build the inequality jacobian
0120             c = [cL; cU];
0121             cJacob = transpose([cJacobL; cJacobU]);
0122         end
0123         
0124         function [statusmsg,statusflag] = GetStatus(obj)
0125             try 
0126                 statusmsg = obj.resultInfo.message;
0127                 if obj.resultInfo.exitflag == 0
0128                     statusflag = 0;
0129                 elseif obj.resultInfo.exitflag > 0
0130                     statusflag = 1;
0131                 else 
0132                     statusflag = -1;
0133                 end
0134             catch
0135                 statusmsg = 'No Last Optimizer Status Available!';
0136                 statusflag = -1;
0137             end
0138         end
0139     end
0140     
0141     methods (Static)    
0142         function available = IsAvailable()
0143             %'fmincon' is a p-code file in the optimization toolbox
0144             available = exist('fmincon') == 6;
0145         end
0146     end
0147 end