diff options
Diffstat (limited to 'macros/fminunc.sci')
| -rw-r--r-- | macros/fminunc.sci | 427 |
1 files changed, 427 insertions, 0 deletions
diff --git a/macros/fminunc.sci b/macros/fminunc.sci new file mode 100644 index 0000000..07fa92a --- /dev/null +++ b/macros/fminunc.sci @@ -0,0 +1,427 @@ +// Copyright (C) 2015 - IIT Bombay - FOSSEE +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Author: R.Vidyadhar & Vignesh Kannan +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in + +function [xopt,fopt,exitflag,output,gradient,hessian] = fminunc (varargin) + // Solves a multi-variable unconstrainted optimization problem + // + // Calling Sequence + // xopt = fminunc(f,x0) + // xopt = fminunc(f,x0,options) + // [xopt,fopt] = fminunc(.....) + // [xopt,fopt,exitflag]= fminunc(.....) + // [xopt,fopt,exitflag,output]= fminunc(.....) + // [xopt,fopt,exitflag,output,gradient]=fminunc(.....) + // [xopt,fopt,exitflag,output,gradient,hessian]=fminunc(.....) + // + // Parameters + // f : a function, representing the objective function of the problem + // x0 : a vector of doubles, containing the starting of variables. + // options: a list, containing the option for user to specify. See below for details. + // xopt : a vector of doubles, the computed solution of the optimization problem. + // fopt : a scalar of double, the function value at x. + // exitflag : a scalar of integer, containing the flag which denotes the reason for termination of algorithm. See below for details. + // output : a structure, containing the information about the optimization. See below for details. + // gradient : a vector of doubles, containing the the gradient of the solution. + // hessian : a matrix of doubles, containing the the hessian of the solution. + // + // Description + // Search the minimum of an unconstrained optimization problem specified by : + // Find the minimum of f(x) such that + // + // <latex> + // \begin{eqnarray} + // &\mbox{min}_{x} + // & f(x)\\ + // \end{eqnarray} + // </latex> + // + // The routine calls Ipopt for solving the Un-constrained Optimization problem, Ipopt is a library written in C++. + // + // The options allows the user to set various parameters of the Optimization problem. + // It should be defined as type "list" and contains the following fields. + // <itemizedlist> + // <listitem>Syntax : options= list("MaxIter", [---], "CpuTime", [---], "Gradient", ---, "Hessian", ---);</listitem> + // <listitem>MaxIter : a Scalar, containing the Maximum Number of Iteration that the solver should take.</listitem> + // <listitem>CpuTime : a Scalar, containing the Maximum amount of CPU Time that the solver should take.</listitem> + // <listitem>Gradient : a function, representing the gradient function of the Objective in Vector Form.</listitem> + // <listitem>Hessian : a function, representing the hessian function of the Objective in Symmetric Matrix Form.</listitem> + // <listitem>Default Values : options = list("MaxIter", [3000], "CpuTime", [600]);</listitem> + // </itemizedlist> + // + // The exitflag allows to know the status of the optimization which is given back by Ipopt. + // <itemizedlist> + // <listitem>exitflag=0 : Optimal Solution Found </listitem> + // <listitem>exitflag=1 : Maximum Number of Iterations Exceeded. Output may not be optimal.</listitem> + // <listitem>exitflag=2 : Maximum CPU Time exceeded. Output may not be optimal.</listitem> + // <listitem>exitflag=3 : Stop at Tiny Step.</listitem> + // <listitem>exitflag=4 : Solved To Acceptable Level.</listitem> + // <listitem>exitflag=5 : Converged to a point of local infeasibility.</listitem> + // </itemizedlist> + // + // For more details on exitflag see the ipopt documentation, go to http://www.coin-or.org/Ipopt/documentation/ + // + // The output data structure contains detailed informations about the optimization process. + // It has type "struct" and contains the following fields. + // <itemizedlist> + // <listitem>output.Iterations: The number of iterations performed during the search</listitem> + // <listitem>output.Cpu_Time: The total cpu-time spend during the search</listitem> + // <listitem>output.Objective_Evaluation: The number of Objective Evaluations performed during the search</listitem> + // <listitem>output.Dual_Infeasibility: The Dual Infeasiblity of the final soution</listitem> + // </itemizedlist> + // + // Examples + // //Find x in R^2 such that it minimizes the Rosenbrock function + // //f = 100*(x2 - x1^2)^2 + (1-x1)^2 + // //Objective function to be minimised + // function y= f(x) + // y= 100*(x(2) - x(1)^2)^2 + (1-x(1))^2; + // endfunction + // //Starting point + // x0=[-1,2]; + // //Gradient of objective function + // function y= fGrad(x) + // y= [-400*x(1)*x(2) + 400*x(1)^3 + 2*x(1)-2, 200*(x(2)-x(1)^2)]; + // endfunction + // //Hessian of Objective Function + // function y= fHess(x) + // y= [1200*x(1)^2- 400*x(2) + 2, -400*x(1);-400*x(1), 200 ]; + // endfunction + // //Options + // options=list("MaxIter", [1500], "CpuTime", [500], "Gradient", fGrad, "Hessian", fHess); + // //Calling Ipopt + // [xopt,fopt,exitflag,output,gradient,hessian]=fminunc(f,x0,options) + // // Press ENTER to continue + // + // Examples + // //Find x in R^2 such that the below function is minimum + // //f = x1^2 + x2^2 + // //Objective function to be minimised + // function y= f(x) + // y= x(1)^2 + x(2)^2; + // endfunction + // //Starting point + // x0=[2,1]; + // //Calling Ipopt + // [xopt,fopt]=fminunc(f,x0) + // // Press ENTER to continue + // + // Examples + // //The below problem is an unbounded problem: + // //Find x in R^2 such that the below function is minimum + // //f = - x1^2 - x2^2 + // //Objective function to be minimised + // function y= f(x) + // y= -x(1)^2 - x(2)^2; + // endfunction + // //Starting point + // x0=[2,1]; + // //Gradient of objective function + // function y= fGrad(x) + // y= [-2*x(1),-2*x(2)]; + // endfunction + // //Hessian of Objective Function + // function y= fHess(x) + // y= [-2,0;0,-2]; + // endfunction + // //Options + // options=list("MaxIter", [1500], "CpuTime", [500], "Gradient", fGrad, "Hessian", fHess); + // //Calling Ipopt + // [xopt,fopt,exitflag,output,gradient,hessian]=fminunc(f,x0,options) + // Authors + // R.Vidyadhar , Vignesh Kannan + + + //To check the number of input and output arguments + [lhs , rhs] = argn(); + + //To check the number of arguments given by the user + if ( rhs<2 | rhs>5 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2 or 5"), "fminunc", rhs); + error(errmsg) + end + + //Storing the 1st and 2nd Input Parameters + fun = varargin(1); + x0 = varargin(2); + + //To check whether the 1st Input argument(f) is a function or not + if (type(f) ~= 13 & type(f) ~= 11) then + errmsg = msprintf(gettext("%s: Expected function for Objective "), "fminunc"); + error(errmsg); + end + + //To check whether the 2nd Input argument(x0) is a vector/scalar + if (type(x0) ~= 1) then + errmsg = msprintf(gettext("%s: Expected Vector/Scalar for Starting Point"), "fminunc"); + error(errmsg); + end + + //To check and convert the 2nd Input argument(x0) to a row vector + if((size(x0,1)~=1) & (size(x0,2)~=1)) then + errmsg = msprintf(gettext("%s: Expected Row Vector or Column Vector for x0 (Initial Value) "), "fminunc", rhs); + error(errmsg); + else + if(size(x0,2)==1) then + x0=x0'; //Converting x0 to a row vector, if it is a column vector + else + x0=x0; //Retaining the same, if it is already a row vector + end + s=size(x0); + end + + + //To check the match between f (1st Parameter) and x0 (2nd Parameter) + if(execstr('init=fun(x0)','errcatch')==21) then + errmsg = msprintf(gettext("%s: Objective function and x0 did not match"), "fminunc"); + error(errmsg); + end + + //Converting the User defined Objective function into Required form (Error Detectable) + function [y,check] = f(x) + if(execstr('y=fun(x)','errcatch')==32 | execstr('y=fun(x)','errcatch')==27) + y=0; + check=1; + else + y=fun(x); + if (isreal(y)==%F) then + y=0; + check=1; + else + check=0; + end + end + endfunction + + //To check whether options has been entered by user + if ( rhs<3 ) then + param = list(); + else + param =varargin(3); //Storing the 3rd Input Parameter in an intermediate list named 'param' + + end + + //If options has been entered, then check its type for 'list' + if (type(param) ~= 15) then + errmsg = msprintf(gettext("%s: 3rd Input parameter should be a list (ie. Options) "), "fminunc"); + error(errmsg); + end + + //If options has been entered, then check whether an even number of entires has been entered + if (modulo(size(param),2)) then + errmsg = msprintf(gettext("%s: Size of parameters should be even"), "fminunc"); + error(errmsg); + end + + //Defining a function to calculate Gradient or Hessian if the respective user entry is OFF + function [y,check]=gradhess(x,t) + if t==1 then //To return Gradient + if(execstr('y=numderivative(fun,x)','errcatch')==10000) + y=0; + check=1; + else + y=numderivative(fun,x); + if (isreal(y)==%F) then + y=0; + check=1; + else + check=0; + end + end + else //To return Hessian + if(execstr('[grad,y]=numderivative(fun,x)','errcatch')==10000) + y=0; + check=1; + else + [grad,y]=numderivative(fun,x); + if (isreal(y)==%F) then + y=0; + check=1; + else + check=0; + end + end + end + endfunction + + //To set default values for options, if user doesn't enter options + options = list("MaxIter", [3000], "CpuTime", [600]); + + //Flags to check whether Gradient is "ON"/"OFF" and Hessian is "ON"/"OFF" + flag1=0; + flag2=0; + fGrad=[]; + fGrad1=[]; + fHess=[]; + fHess1=[]; + + //To check the user entry for options and store it + for i = 1:(size(param))/2 + select param(2*i-1) + case "MaxIter" then + options(2*i) = param(2*i); //Setting the maximum number of iterations as per user entry + case "CpuTime" then + options(2*i) = param(2*i); //Setting the maximum CPU time as per user entry + case "Gradient" then + flag1 = 1; + fGrad = param(2*i); + case "Hessian" then + flag2 = 1; + fHess = param(2*i); + else + errmsg = msprintf(gettext("%s: Unrecognized parameter name %s."), "fminbnd", param(2*i-1)); + error(errmsg) + end + end + + + //To check for correct input of Gradient and Hessian functions from the user + if (flag1==1) then + if (type(fGrad) ~= 13 & type(fGrad) ~= 11) then + errmsg = msprintf(gettext("%s: Expected function for Gradient of Objective"), "fminunc"); + error(errmsg); + end + + if(execstr('samplefGrad=fGrad(x0)','errcatch')==21) + errmsg = msprintf(gettext("%s: Gradient function of Objective and x0 did not match"), "fminunc"); + error(errmsg); + end + + samplefGrad=fGrad(x0); + + if (size(samplefGrad,1)==s(2) & size(samplefGrad,2)==1) then + elseif (size(samplefGrad,1)==1 & size(samplefGrad,2)==s(2)) then + elseif (size(samplefGrad,1)~=1 & size(samplefGrad,2)~=1) then + errmsg = msprintf(gettext("%s: Wrong Input for Objective Gradient function(3rd Parameter)---->Row Vector function is Expected"), "fminunc"); + error(errmsg); + end + + function [y,check] = fGrad1(x) + if(execstr('y=fGrad(x)','errcatch')==32 | execstr('y=fGrad(x)','errcatch')==27) + y = 0; + check=1; + else + y=fGrad(x); + if (isreal(y)==%F) then + y = 0; + check=1; + else + check=0; + end + end + endfunction + end + if (flag2==1) then + if (type(fHess) ~= 13 & type(fHess) ~= 11) then + errmsg = msprintf(gettext("%s: Expected function for Hessian of Objective"), "fminunc"); + error(errmsg); + end + + if(execstr('samplefHess=fHess(x0)','errcatch')==21) + errmsg = msprintf(gettext("%s: Hessian function of Objective and x0 did not match"), "fminunc"); + error(errmsg); + end + + samplefHess=fHess(x0); + + if(size(samplefHess,1)~=s(2) | size(samplefHess,2)~=s(2)) then + errmsg = msprintf(gettext("%s: Wrong Input for Objective Hessian function(3rd Parameter)---->Symmetric Matrix function is Expected "), "fminunc"); + error(errmsg); + end + + function [y,check] = fHess1(x) + if(execstr('y=fHess(x)','errcatch')==32 | execstr('y=fHess(x)','errcatch')==27) + y = 0; + check=1; + else + y=fHess(x); + if (isreal(y)==%F) then + y = 0; + check=1; + else + check=0; + end + end + endfunction + end + + //Calling the Ipopt function for solving the above problem + [xopt,fopt,status,iter,cpu,obj_eval,dual,gradient, hessian1] = solveminuncp(f,gradhess,flag1,fGrad1,flag2,fHess1,x0,options); + + //Calculating the values for output + xopt = xopt'; + exitflag = status; + output = struct("Iterations", [],"Cpu_Time",[],"Objective_Evaluation",[],"Dual_Infeasibility",[]); + output.Iterations = iter; + output.Cpu_Time = cpu; + output.Objective_Evaluation = obj_eval; + output.Dual_Infeasibility = dual; + + //Converting hessian of order (1 x (numberOfVariables)^2) received from Ipopt to order (numberOfVariables x numberOfVariables) + s=size(gradient) + for i =1:s(2) + for j =1:s(2) + hessian(i,j)= hessian1(j+((i-1)*s(2))) + end + end + + + //In the cases of the problem not being solved, return NULL to the output matrices + if( status~=0 & status~=1 & status~=2 & status~=3 & status~=4 & status~=7 ) then + xopt=[] + fopt=[] + output = struct("Iterations", [],"Cpu_Time",[]); + output.Iterations = iter; + output.Cpu_Time = cpu; + gradient=[] + hessian=[] + end + + + //To print output message + select status + + case 0 then + printf("\nOptimal Solution Found.\n"); + case 1 then + printf("\nMaximum Number of Iterations Exceeded. Output may not be optimal.\n"); + case 2 then + printf("\nMaximum CPU Time exceeded. Output may not be optimal.\n"); + case 3 then + printf("\nStop at Tiny Step\n"); + case 4 then + printf("\nSolved To Acceptable Level\n"); + case 5 then + printf("\nConverged to a point of local infeasibility.\n"); + case 6 then + printf("\nStopping optimization at current point as requested by user.\n"); + case 7 then + printf("\nFeasible point for square problem found.\n"); + case 8 then + printf("\nIterates diverging; problem might be unbounded.\n"); + case 9 then + printf("\nRestoration Failed!\n"); + case 10 then + printf("\nError in step computation (regularization becomes too large?)!\n"); + case 12 then + printf("\nProblem has too few degrees of freedom.\n"); + case 13 then + printf("\nInvalid option thrown back by Ipopt\n"); + case 14 then + printf("\nNot enough memory.\n"); + case 15 then + printf("\nINTERNAL ERROR: Unknown SolverReturn value - Notify Ipopt Authors.\n"); + else + printf("\nInvalid status returned. Notify the Toolbox authors\n"); + break; + end + + +endfunction |