diff options
Diffstat (limited to 'macros')
38 files changed, 3745 insertions, 2 deletions
diff --git a/macros/fgoalattain.bin b/macros/fgoalattain.bin Binary files differnew file mode 100644 index 0000000..1476d8e --- /dev/null +++ b/macros/fgoalattain.bin diff --git a/macros/fgoalattain.sci b/macros/fgoalattain.sci new file mode 100644 index 0000000..5f519b4 --- /dev/null +++ b/macros/fgoalattain.sci @@ -0,0 +1,514 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2015 - IIT Bombay - FOSSEE +// +// Authors: Prajwala TM,Sheetal Shalini +// Organization: FOSSEE, IIT Bombay +// Email: prajwala.tm@gmail.com,sheetalsh456@gmail.com +// 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 + +function [x,fval,attainfactor,exitflag,output,lambda] = fgoalattain(varargin) + // Solves a multiobjective goal attainment problem + // + // Calling Sequence + // x = fgoalattain(fun,x0,goal,weight) + // x = fgoalattain(fun,x0,goal,weight,A,b) + // x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq) + // x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub) + // x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon) + // x = fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon,options) + // [x,fval] = fgoalattain(...) + // [x,fval,attainfactor] = fgoalattain(...) + // [x,fval,attainfactor,exitflag] = fgoalattain(...) + // [x,fval,attainfactor,exitflag,output] = fgoalattain(...) + // [x,fval,attainfactor,exitflag,output,lambda] = fgoalattain(...) + // + // Parameters + // fun: a function that accepts a vector x and returns a vector F + // x0: a nx1 or 1xn matrix of double, where n is the number of variables. + // The initial guess for the optimization algorithm + // A: a nil x n matrix of double, where n is the number of variables and + // nil is the number of linear inequalities. + // + // b: a nil x 1 matrix of double, where nil is the number of linear + // inequalities + // Aeq: a nel x n matrix of double, where n is the number of variables + // and nel is the number of linear equalities. + // beq: a nel x 1 matrix of double, where nel is the number of linear + // equalities + // lb: a nx1 or 1xn matrix of double, where n is the number of variables. + // The lower bound for x. If lb==[], then the lower bound is + // automatically set to -inf + // ub: a nx1 or 1xn matrix of double, where n is the number of variables. + // The upper bound for x. If ub==[], then the upper bound is + // automatically set to +inf + // nonlcon: a function, the nonlinear constraints + // options : a list, containing the option for user to specify. See below for details. + // x: a nx1 matrix of double, the computed solution of the optimization problem + // fval: a vector of double, the value of functions at x + // attainfactor: The amount of over- or underachievement of the goals,γ at the solution. + // exitflag: a 1x1 matrix of floating point integers, the exit status + // output: a struct, the details of the optimization process + // lambda: a struct, the Lagrange multipliers at optimum + // + // Description + // fgoalattain solves the goal attainment problem, which is one formulation for minimizing a multiobjective optimization problem. + // Finds the minimum of a problem specified by: + // Minimise Y such that + // + //<latex> + //\begin{eqnarray} + //\mbox{min}_{x,\gamma} & f(x)-weight \ast \gamma \leq goal \\ + //\mbox{subject to} & c(x) \leq 0 \\ + // & c_{eq}(x) = 0 \\ + // & Ax \leq b \\ + // & A_{eq} x = b_{eq} \\ + // & lb \leq x \leq ub + //\end{eqnarray} + //</latex> + // + // The solver makes use of fmincon to find the minimum. + // + // The fgoalattain finds out the maximum value of Y for the objectives evaluated at the starting point and + // adds that as another variable to the vector x + // This is passed to the fmincon function to get the optimised value of Y + // Hence, the algorithm used mainly is "ipopt" to obtain the optimum solution + // The relations between f(x), Y, weights and goals are added as additional non-linear inequality constraints + // + // 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", [---], "GradObj", ---, "GradCon", ---);</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>GradObj : a function, representing the gradient function of the Objective in Vector Form.</listitem> + // <listitem>GradCon : a function, representing the gradient of the Non-Linear Constraints (both Equality and Inequality) of the problem. It is declared in such a way that gradient of non-linear inequality constraints are defined first as a separate Matrix (cg of size m2 X n or as an empty), followed by gradient of non-linear equality constraints as a separate Matrix (ceqg of size m2 X n or as an empty) where m2 & m3 are number of non-linear inequality and equality constraints respectively.</listitem> + // <listitem>Default Values : options = list("MaxIter", [3000], "CpuTime", [600]);</listitem> + // </itemizedlist> + // + // By default, the gradient options for fminimax are turned off and and fmincon does the gradient opproximation of minmaxObjfun. In case the GradObj option is off and GradConstr option is on, fminimax approximates minmaxObjfun gradient using numderivative toolbox. + // + // If we can provide exact gradients, we should do so since it improves the convergence speed of the optimization algorithm. + // + // Furthermore, we must enable the "GradObj" option with the statement : + // <programlisting> + // minimaxOptions = list("GradObj",fGrad); + // </programlisting> + // This will let fminimax know that the exact gradient of the objective function is known, so that it can change the calling sequence to the objective function. Note that, fGrad should be mentioned in the form of N x n where n is the number of variables, N is the number of functions in objective function. + // + // The constraint function must have header : + // <programlisting> + // [c, ceq] = confun(x) + // </programlisting> + // where x is a n x 1 matrix of dominmaxUbles, c is a 1 x nni matrix of doubles and ceq is a 1 x nne matrix of doubles (nni : number of nonlinear inequality constraints, nne : number of nonlinear equality constraints). + // On input, the variable x contains the current point and, on output, the variable c must contain the nonlinear inequality constraints and ceq must contain the nonlinear equality constraints. + // + // By default, the gradient options for fminimax are turned off and and fmincon does the gradient opproximation of confun. In case the GradObj option is on and GradCons option is off, fminimax approximates confun gradient using numderivative toolbox. + // + // If we can provide exact gradients, we should do so since it improves the convergence speed of the optimization algorithm. + // + // Furthermore, we must enable the "GradCon" option with the statement : + // <programlisting> + // minimaxOptions = list("GradCon",confunGrad); + // </programlisting> + // This will let fminimax know that the exact gradient of the objective function is known, so that it can change the calling sequence to the objective function. + // + // The constraint derivative function must have header : + // <programlisting> + // [dc,dceq] = confungrad(x) + // </programlisting> + // where dc is a nni x n matrix of doubles and dceq is a nne x n matrix of doubles. + // + // 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 amount of 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> + // + // The lambda data structure contains the Lagrange multipliers at the end + // of optimization. In the current version the values are returned only when the the solution is optimal. + // It has type "struct" and contains the following fields. + // <itemizedlist> + // <listitem>lambda.lower: The Lagrange multipliers for the lower bound constraints.</listitem> + // <listitem>lambda.upper: The Lagrange multipliers for the upper bound constraints.</listitem> + // <listitem>lambda.eqlin: The Lagrange multipliers for the linear equality constraints.</listitem> + // <listitem>lambda.ineqlin: The Lagrange multipliers for the linear inequality constraints.</listitem> + // <listitem>lambda.eqnonlin: The Lagrange multipliers for the non-linear equality constraints.</listitem> + // <listitem>lambda.ineqnonlin: The Lagrange multipliers for the non-linear inequality constraints.</listitem> + // </itemizedlist> + // + // Examples + // function f1 = fun(x) + // f1(1)=2*x(1)*x(1)+x(2)*x(2)-48*x(1)-40*x(2)+304 + // f1(2)=-x(1)*x(1)-3*x(2)*x(2) + // f1(3)=x(1)+3*x(2)-18 + // f1(4)=-x(1)-x(2) + // f1(5)=x(1)+x(2)-8 + // endfunction + // x0=[-1,1]; + // + // goal=[-5,-3,-2,-1,-4]; + // weight=abs(goal) + // //xopt = [-0.0000011 -63.999998 -2.0000002 -8 3.485D-08] + // //fval = [4 3.99] + // + // //Run fgoalattain + // [xopt,fval,attainfactor,exitflag,output,lambda]=fgoalattain(fun,x0,goal,weight) + // + // Authors + // Prajwala TM, Sheetal Shalini , 2015 + + // Check number of input and output arguments + [gattainLhs,gattainRhs] = argn() + fgoalattainCheckrhs("fgoalattain", gattainRhs, [4 6 8 10 11 12]) + fgoalattainChecklhs("fgoalattain", gattainLhs, 1:6) + + // initialisation of fun + gattainObjfun = varargin(1) + fgoalattainChecktype("fgoalattain", gattainObjfun, "fun", 1, "function") + + // initialisation of x0 + gattainStartpoint = varargin(2) + fgoalattainChecktype("fgoalattain", gattainStartpoint, "x0", 2, "constant") + + gattainNumvar = size(gattainStartpoint,"*") + fgoalattainCheckvector("fgoalattain", gattainStartpoint, "x0", 2, gattainNumvar) + gattainStartpoint = gattainStartpoint(:) + + // initialisation of goal + goal=varargin(3) + fgoalattainChecktype("fgoalattain",goal,"goal",3,"constant") + + // initialisation of weight + weight=varargin(4) + fgoalattainChecktype("fgoalattain",weight,"weight",4,"constant") + + //initialisation of A and b + if(gattainRhs < 5) then + gattainA = [] + gattainB = [] + else + gattainA = varargin(5) + gattainB = varargin(6) + end + + fgoalattainChecktype("fgoalattain", gattainA, "A", 5, "constant") + fgoalattainChecktype("fgoalattain", gattainB, "b", 6, "constant") + + gattainNumrowA = size(gattainA,"r") + if(gattainA <> []) then + fgoalattainCheckdims("fgoalattain", gattainA, "A", 5, [gattainNumrowA gattainNumvar]) + fgoalattainCheckvector("fgoalattain", gattainB, "b", 6, gattainNumrowA) + gattainB = gattainB(:) + end + + //initialisation of Aeq and beq + if(gattainRhs < 7) then + gattainAeq = [] + gattainBeq = [] + else + gattainAeq = varargin(7) + gattainBeq = varargin(8) + end + + fgoalattainChecktype("fgoalattain", gattainAeq, "Aeq", 7, "constant") + fgoalattainChecktype("fgoalattain", gattainBeq, "beq", 8, "constant") + + gattainNumrowAeq = size(gattainAeq,"r") + if(gattainAeq <> []) then + fgoalattainCheckdims("fgoalattain", gattainAeq, "Aeq", 7, [gattainNumrowAeq gattainNumvar]) + fgoalattainCheckvector("fgoalattain", gattainBeq, "beq", 8, gattainNumrowAeq) + gattainBeq = gattainBeq(:) + end + + // initialisation of lb and ub + if(gattainRhs < 9) then + gattainLb = [] + gattainUb = [] + else + gattainLb = varargin(9) + gattainUb = varargin(10) + end + + fgoalattainChecktype("fgoalattain", gattainLb, "lb", 9, "constant") + fgoalattainChecktype("fgoalattain", gattainUb, "ub", 10, "constant") + + // Check dimensions of lb and ub + if(gattainLb <> []) then + fgoalattainCheckvector("fgoalattain", gattainLb, "lb", 9, gattainNumvar) + gattainLb = gattainLb(:) + end + + if(gattainUb <> []) then + fgoalattainCheckvector("fgoalattain", gattainUb, "ub", 10, gattainNumvar) + gattainUb = gattainUb(:) + end + + // Initialisation of nonlcon + function [c,ceq] = constr(z) + c = [] + ceq = [] + endfunction + + if(gattainRhs < 11) then + gattainNonlinfun = constr + else + gattainNonlinfun = varargin(11) + end + + fgoalattainChecktype("fgoalattain", gattainNonlinfun, "nonlcon", 11, "function") + + // initialisation of default options + if(gattainRhs < 12) then + gattainUserOptions = list() + else + gattainUserOptions = varargin(12) + end + + //If gattainOptions is entered then checking its type for 'list' + if (type(gattainUserOptions) ~= 15) then + errmsg = msprintf(gettext("%s: gattainOptions (10th parameter) should be a list"), "fgoalattain"); + error(errmsg); + end + + //If minimaxOptions is entered then checking whether even number of entires are entered + if (modulo(size(gattainUserOptions),2)) then + errmsg = msprintf(gettext("%s: Size of gattainOptions (list) should be even"), "fgoalattain"); + error(errmsg); + end + + //Flags to check whether Gradient is "ON"/"OFF" & Hessian is "ON"/"OFF" + flag1=0; + flag2=0; + fgaMaxIter = 3000; + fgaCPU = 600; + gattainFGrad=[]; + gattainCGrad=[]; + + //To check the User Entry for Options and storing it + for i = 1:(size(gattainUserOptions))/2 + select gattainUserOptions(2*i-1) + case "MaxIter" then + fgaMaxIter = gattainUserOptions(2*i); //Setting the Maximum Iteration as per user entry + + case "CpuTime" then + fgaCPU = gattainUserOptions(2*i); //Setting the Maximum CPU Time as per user entry + + case "GradObj" then + flag1=1; + gattainFGrad = gattainUserOptions(2*i); + + case "GradCon" then + flag2=1; + gattainCGrad = gattainUserOptions(2*i); + + else + errmsg = msprintf(gettext("%s: Unrecognized gattainUserOptionseter name ''%s''."), "fminimax", gattainUserOptions(2*i-1)); + error(errmsg) + end + end + + // Checking if minmaxFGrad and minmaxCGrad are functions + if (flag1==1) then + if (type(gattainFGrad) ~= 11 & type(gattainFGrad) ~= 13) then + errmsg = msprintf(gettext("%s: Expected function for Gradient of Objective"), "fminimax"); + error(errmsg); + end + end + + if (flag2==1) then + if (type(gattainCGrad) ~= 11 & type(gattainCGrad) ~= 13) then + errmsg = msprintf(gettext("%s: Expected function for Gradient of Nonlinfun"), "fminimax"); + error(errmsg); + end + end + + gattainObjfunval = gattainObjfun(gattainStartpoint) + gattainObjfunval=gattainObjfunval(:) + goal=goal(:) + weight=weight(:) + gaVal=[] + + // appending the gamma value as another variable + for i=1:size(gattainObjfunval,'r') + if(weight(i)<>0) then + gaVal(i)=((gattainObjfunval(i)-goal(i))/weight(i)) + end + end + + gattainStartpoint(gattainNumvar+1)=max(gaVal) + + if(gattainA <> []) then + gattainA = [gattainA, zeros(gattainNumrowA,1)] + end + if(gattainAeq <> []) then + gattainAeq = [gattainAeq, zeros(gattainNumrowAeq,1)] + end + if(gattainLb <> []) then + gattainLb(gattainNumvar+1) = -%inf + end + if(gattainUb <> []) then + gattainUb(gattainNumvar+1) = +%inf + end + + // function handle defining the additional inequalities + function temp = gattainAddIneq(z) + gaVar = gattainObjfun(z) + gattainAddIneqWithWt = [] + gattainAddIneqWitoutWt = [] + for i = 1:size(gaVar,'r') + if(weight(i) <> 0) then + gattainAddIneqWithWt = [gattainAddIneqWithWt; ( (gaVar(i)-goal(i))/weight(i) )] + else + gattainAddIneqWitoutWt = [gattainAddIneqWitoutWt; gaVar(i)-goal(i)] + end + end + temp = [gattainAddIneqWithWt - z(gattainNumvar+1); gattainAddIneqWitoutWt] + endfunction + + // function handle defining new objective function + function newfunc = newObjfun(z) + newfunc = z(gattainNumvar+1) + endfunction + + // function handle defining add_ineq derivative using numderivative + function func = gattainIneqDer(z) + func = numderivative(gattainAddIneq,z) + endfunction + + // function handle defining nonlcon derivative using numderivative + function [dc,dceq] = gattainNonlinDer(z) + + // function handle extracting c and ceq components from nonlcon + function foo = gattainC(z) + [foo,tmp1] = gattainNonlinfun(z) + foo = foo' + endfunction + + function foo = gattainCEQ(z) + [tmp1,foo] = gattainNonlinfun(z) + foo = foo' + endfunction + + dc = numderivative(gattainC,z) + dceq = numderivative(gattainCEQ,z) + endfunction + + // function handle defining new nonlcon function + function [nc,nceq] = newNonlinfun(z) + [nc,nceq] = gattainNonlinfun(z) + tmp = [gattainAddIneq(z)]' + nc = [nc, tmp] + endfunction + + function [dnc,dnceq] = newCGrad(z) + // check if "GradCon" option is turned on + // if "GradCon" is turned on, use it + if(flag2 == 1) then + [dnc,dnceq] = gattainCGrad(z) + dnc = [dnc, zeros(size(dnc,'r'),1)] + dnceq = [dnceq, zeros(size(dnceq,'r'),1)] + // else, calculate it using finite differences + else + [dnc,dnceq] = gattainNonlinDer(z) + end + + // check if "GradObj" option is turned on + // if "GradObj" is turned on, use it + if(flag1 == 1) then + derObjfun = gattainFGrad(z) + tmp1 = [] + tmp2 = [] + for i = 1:size(gattainObjfun(gattainStartpoint),'r') + if weight(i) <> 0 then + gaVal = [derObjfun(i,:)/weight(i) , -1] + tmp1 = [ tmp1; gaVal ] + else + gaVal = [derObjfun(i,:) , 0] + tmp2 = [ tmp2; gaVal ] + end + end + + dnc = [ dnc; tmp1; tmp2 ] + // else, calculate it using finite differences + else + deraddineq = gattainIneqDer(z) + dnc = [dnc; deraddineq] + end + endfunction + + + // disp(gattainStartpoint) + // disp(gattainObjfun(gattainStartpoint)) + // disp(newObjfun(gattainStartpoint)) + // disp(goal) + // disp(weight) + // + // disp(gattainA) + // disp(gattainB) + // disp(gattainAeq) + // disp(gattainBeq) + // disp(gattainLb) + // disp(gattainUb) + // + // [f,g] = gattainNonlinfun(gattainStartpoint) + // disp(f) + // disp(g) + // [f,g] = newNonlinfun(gattainStartpoint) + // disp(f) + // disp(g) + // + // // function foo = CgattainC(z) + // // [foo,tmp1] = newNonlinfun(z) + // // endfunction + // // + // // function foo = CgattainCEQ(z) + // // [tmp1,foo] = newNonlinfun(z) + // // endfunction + // // + // // // function handle defining gattainNonlinfun derivative using numderivative + // // function [dc,dceq] = derrNonlinApp(z) + // // dc = numderivative(CgattainC,z) + // // dceq = numderivative(CgattainCEQ,z) + // // endfunction + // + // // [f,g] = derrNonlinApp(gattainStartpoint) + // // disp(f) + // // disp(g) + // + // [f,g] = newCGrad(gattainStartpoint) + // disp(f) + // disp(g) + + //to be passed as options to fmincon + if (flag1 == 1 | flag2 == 1) then + gattainPassOptions = list("MaxIter", fgaMaxIter, "CpuTime", fgaCPU, "GradCon", newCGrad) + [x,attainfactor,exitflag,output,lambda] = fmincon(newObjfun,gattainStartpoint,gattainA,gattainB,gattainAeq,gattainBeq,gattainLb,gattainUb,newNonlinfun,gattainPassOptions) + x= x(1:gattainNumvar) + fval = gattainObjfun(x) + else + gattainPassOptions = list("MaxIter", fgaMaxIter, "CpuTime", fgaCPU) + [x,attainfactor,exitflag,output,lambda] = fmincon(newObjfun,gattainStartpoint,gattainA,gattainB,gattainAeq,gattainBeq,gattainLb,gattainUb,newNonlinfun,gattainPassOptions) + x= x(1:gattainNumvar) + fval = gattainObjfun(x) + end + +endfunction diff --git a/macros/fgoalattainFunctions.bin b/macros/fgoalattainFunctions.bin Binary files differnew file mode 100644 index 0000000..eb72f4d --- /dev/null +++ b/macros/fgoalattainFunctions.bin diff --git a/macros/fgoalattainFunctions.sci b/macros/fgoalattainFunctions.sci new file mode 100644 index 0000000..8745fa6 --- /dev/null +++ b/macros/fgoalattainFunctions.sci @@ -0,0 +1,76 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2015 - IIT Bombay - FOSSEE +// +// Authors: Prajwala TM,Sheetal Shalini +// Organization: FOSSEE, IIT Bombay +// Email: prajwala.tm@gmail.com,sheetalsh456@gmail.com +// 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 + +// Define subsidiary functions used by fgoalattain + +// to check if the number of input arguments is equal to the numbers mentioned in the set +function errmsg = fgoalattainCheckrhs ( funname , rhs , rhsset ) + errmsg = [] + if ( and(rhs <> rhsset) ) then + rhsstr = strcat(string(rhsset)," ") + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while the number of expected input arguments should be in the set [%s]."), funname , rhs , rhsstr ); + error(errmsg) + end +endfunction + +// to check if the number of output arguments is equal to the numbers mentioned in the set, which is any number between 1 to 6 (both inclusive) +function errmsg = fgoalattainChecklhs ( funname , lhs , lhsset ) +errmsg = [] + if ( and ( lhs <> lhsset ) ) then + lhsstr = strcat(string(lhsset)," ") + errmsg = msprintf(gettext("%s: Unexpected number of output arguments : %d provided while the expected number of output arguments should be in the set [%s]."), funname , lhs , lhsstr ); + error(errmsg) + end +endfunction + +// Generates an error if the given variable is not of expected type. +function errmsg = fgoalattainChecktype ( funname , var , varname , ivar , expectedtype ) +errmsg = [] + if ( and ( typeof ( var ) <> expectedtype ) ) then + strexp = """" + strcat(expectedtype,""" or """) + """" + errmsg = msprintf(gettext("%s: Expected type [%s] for input argument %s at input #%d, but got ""%s"" instead."),funname, strexp, varname , ivar , typeof(var) ); + error(errmsg); + end +endfunction + +// Generates an error if the variable is not a vector. +function errmsg = fgoalattainCheckvector ( funname , var , varname , ivar , nbval ) +errmsg = [] + nrows = size(var,"r") + ncols = size(var,"c") + if ( nrows <> 1 & ncols <> 1 ) then + strcomp = strcat(string(size(var))," ") + errmsg = msprintf(gettext("%s: Expected a vector matrix for input argument %s at input #%d, but got [%s] instead."), funname, varname , ivar , strcomp ); + error(errmsg) + end + if ( ( nrows == 1 & ncols <> nbval ) | ( ncols == 1 & nrows <> nbval ) ) then + strcomp = strcat(string(size(var))," ") + errmsg = msprintf(gettext("%s: Expected %d entries for input argument %s at input #%d, but current dimensions are [%s] instead."), funname, nbval , varname , ivar , strcomp ); + error(errmsg) + end +endfunction + +// Generates an error if the variable has not the required size. +function errmsg = fgoalattainCheckdims ( funname , var , varname , ivar , matdims ) +[lhs,rhs]=argn() + fgoalattainCheckrhs ( "fgoalattainCheckdims" , rhs , 5 ) + fgoalattainChecklhs ( "fgoalattainCheckdims" , lhs , [0 1] ) + errmsg = [] + if ( or ( size(var) <> matdims ) ) then + strexp = strcat(string(matdims)," ") + strcomp = strcat(string(size(var))," ") + errmsg = msprintf(gettext("%s: Expected size [%s] for input argument %s at input #%d, but got [%s] instead."), funname, strexp, varname , ivar , strcomp ); + error(errmsg) + end +endfunction + + diff --git a/macros/fminbnd.bin b/macros/fminbnd.bin Binary files differnew file mode 100644 index 0000000..76db1fb --- /dev/null +++ b/macros/fminbnd.bin diff --git a/macros/fminbnd.sci b/macros/fminbnd.sci new file mode 100644 index 0000000..2ebaa8c --- /dev/null +++ b/macros/fminbnd.sci @@ -0,0 +1,356 @@ +// 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,lambda] = fminbnd (varargin) + // Solves a multi-variable optimization problem on a bounded interval + // + // Calling Sequence + // xopt = fminbnd(f,x1,x2) + // xopt = fminbnd(f,x1,x2,options) + // [xopt,fopt] = fminbnd(.....) + // [xopt,fopt,exitflag]= fminbnd(.....) + // [xopt,fopt,exitflag,output]=fminbnd(.....) + // [xopt,fopt,exitflag,output,lambda]=fminbnd(.....) + // + // Parameters + // f : a function, representing the objective function of the problem + // x1 : a vector, containing the lower bound of the variables of size (1 X n) or (n X 1) where 'n' is the number of Variables, where n is number of Variables + // x2 : a vector, containing the upper bound of the variables of size (1 X n) or (n X 1) or (0 X 0) where 'n' is the number of Variables. If x2 is empty it means upper bound is +infinity + // options : a list, containing the option for user to specify. See below for details. + // xopt : a vector of doubles, containing the the computed solution of the optimization problem. + // fopt : a scalar of double, containing the 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. + // lambda : a structure, containing the Lagrange multipliers of lower bound and upper bound at the optimized point. See below for details. + // + // Description + // Search the minimum of a multi-variable function on bounded interval specified by : + // Find the minimum of f(x) such that + // + // <latex> + // \begin{eqnarray} + // &\mbox{min}_{x} + // & f(x)\\ + // & \text{subject to} & x1 \ < x \ < x2 \\ + // \end{eqnarray} + // </latex> + // + // The routine calls Ipopt for solving the Bounded 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", [---], TolX, [----]);</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>TolX : a Scalar, containing the Tolerance value that the solver should take.</listitem> + // <listitem>Default Values : options = list("MaxIter", [3000], "CpuTime", [600], TolX, [1e-4]);</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> + // + // The lambda data structure contains the Lagrange multipliers at the end + // of optimization. In the current version the values are returned only when the the solution is optimal. + // It has type "struct" and contains the following fields. + // <itemizedlist> + // <listitem>lambda.lower: The Lagrange multipliers for the lower bound constraints.</listitem> + // <listitem>lambda.upper: The Lagrange multipliers for the upper bound constraints.</listitem> + // </itemizedlist> + // + // Examples + // //Find x in R^6 such that it minimizes: + // //f(x)= sin(x1) + sin(x2) + sin(x3) + sin(x4) + sin(x5) + sin(x6) + // //-2 <= x1,x2,x3,x4,x5,x6 <= 2 + // //Objective function to be minimised + // function y=f(x) + // y=0 + // for i =1:6 + // y=y+sin(x(i)); + // end + // endfunction + // //Variable bounds + // x1 = [-2, -2, -2, -2, -2, -2]; + // x2 = [2, 2, 2, 2, 2, 2]; + // //Options + // options=list("MaxIter",[1500],"CpuTime", [100],"TolX",[1e-6]) + // //Calling Ipopt + // [x,fval] =fminbnd(f, x1, x2, options) + // // Press ENTER to continue + // + // Examples + // //Find x in R such that it minimizes: + // //f(x)= 1/x^2 + // //0 <= x <= 1000 + // //Objective function to be minimised + // function y=f(x) + // y=1/x^2 + // endfunction + // //Variable bounds + // x1 = [0]; + // x2 = [1000]; + // //Calling Ipopt + // [x,fval,exitflag,output,lambda] =fminbnd(f, x1, x2) + // // Press ENTER to continue + // + // Examples + // //The below problem is an unbounded problem: + // //Find x in R^2 such that it minimizes: + // //f(x)= -[(x1-1)^2 + (x2-1)^2] + // //-inf <= x1,x2 <= inf + // //Objective function to be minimised + // function y=f(x) + // y=-((x(1)-1)^2+(x(2)-1)^2); + // endfunction + // //Variable bounds + // x1 = [-%inf , -%inf]; + // x2 = []; + // //Options + // options=list("MaxIter",[1500],"CpuTime", [100],"TolX",[1e-6]) + // //Calling Ipopt + // [x,fval,exitflag,output,lambda] =fminbnd(f, x1, x2, 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<3 | rhs>4 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 3 or 4"), "fminbnd", rhs); + error(errmsg) + end + + //Storing the 1st and 2nd Input Parameters + fun = varargin(1); + x1 = varargin(2); + x2 = varargin(3); + + //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=zeros(1,1); + check=1; + else + if (isreal(y)==%F) then + y=zeros(1,1); + check=1; + else + y=fun(x); + check=0; + end + end + endfunction + + //To check whether the 2nd Input argument (x1) is a vector/scalar + if (type(x1) ~= 1) then + errmsg = msprintf(gettext("%s: Expected Vector/Scalar for Lower Bound Vector (2nd Parameter)"), "fminbnd"); + error(errmsg); + end + + if (size(x1,2)==0) then + errmsg = msprintf(gettext("%s: Lower Bound (2nd Parameter) cannot be empty"), "fminbnd"); + error(errmsg); + end + + if (size(x1,1)~=1) & (size(x1,2)~=1) then + errmsg = msprintf(gettext("%s: Lower Bound (2nd Parameter) should be a vector"), "fminbnd"); + error(errmsg); + elseif (size(x1,2)==1) then + x1=x1; + elseif (size(x1,1)==1) then + x1=x1'; + end + s=size(x1); + + //To check whether the 3rd Input argument (x2) is a vector/scalar + if (type(x2) ~= 1) then + errmsg = msprintf(gettext("%s: Expected Vector/Scalar for Upper Bound Vector (3rd Parameter)"), "fminbnd"); + error(errmsg); + end + + //To check for the correct size and data of x2 (3rd paramter) and convert it to a column vector as required by Ipopt + if (size(x2,2)==0) then + x2 = repmat(%inf,s(1),1); + end + + + if (size(x2,1)~=1) & (size(x2,2)~=1) then + errmsg = msprintf(gettext("%s: Upper Bound (3rd Parameter) should be a vector"), "fminbnd"); + error(errmsg); + elseif(size(x2,1)~=s(1) & size(x2,2)==1) then + errmsg = msprintf(gettext("%s: Upper Bound and Lower Bound are not matching"), "fminbnd"); + error(errmsg); + elseif(size(x2,1)==s(1) & size(x2,2)==1) then + x2=x2; + elseif(size(x2,1)==1 & size(x2,2)~=s(1)) then + errmsg = msprintf(gettext("%s: Upper Bound and Lower Bound are not matching"), "fminbnd"); + error(errmsg); + elseif(size(x2,1)==1 & size(x2,2)==s(1)) then + x2=x2'; + end + + //To check the contents of x1 and x2 (2nd & 3rd Parameter) + + for i = 1:s(2) + if (x1(i) == %inf) then + errmsg = msprintf(gettext("%s: Value of Lower Bound can not be infinity"), "fminbnd"); + error(errmsg); + end + if (x2(i) == -%inf) then + errmsg = msprintf(gettext("%s: Value of Upper Bound can not be negative infinity"), "fminbnd"); + error(errmsg); + end + if(x2(i)-x1(i)<=1e-6) then + errmsg = msprintf(gettext("%s: Difference between Upper Bound and Lower bound should be atleast > 10^6 for variable number= %d "), "fminbnd", i); + error(errmsg) + end + end + + //To check, whether options has been entered by the user + if ( rhs<4 | size(varargin(4)) ==0 ) then + param = list(); + else + param =varargin(4); //Storing the 3rd Input Parameter in an intermediate list named 'param' + end + + options = list("MaxIter",[3000],"CpuTime",[600],"TolX", [1e-4]); + + //To check the user entry for options and storing it + for i = 1:(size(param))/2 + select param(2*i-1) + case "MaxIter" then + options(2*i) = param(2*i); + case "CpuTime" then + options(2*i) = param(2*i); + case "TolX" then + options(2*i) = param(2*i); + else + errmsg = msprintf(gettext("%s: Unrecognized parameter name %s."), "fminbnd", param(2*i-1)); + error(errmsg) + end + end + + + //Defining a function to calculate Gradient or Hessian in an error deductable form. + function [y,check]=gradhess(x,t) + if t==1 then + if(execstr('y=numderivative(fun,x)','errcatch')==10000) + y=zeros(s(1),1); + check=1; + else + y=numderivative(fun,x); + if (isreal(y)==%F) then + y=zeros(s(1),1); + check=1; + else + check=0; + end + end + else + if(execstr('[grad,y]=numderivative(fun,x)','errcatch')==10000) + y=zeros(s(1),1); + check=1; + else + [grad,y]=numderivative(fun,x); + if (isreal(y)==%F) then + y=zeros(s(1),s(1)); + check=1; + else + check=0; + end + end + end + endfunction + + //Calling the Ipopt function for solving the above problem + [xopt,fopt,status,iter,cpu,obj_eval,dual,zl,zu] = solveminbndp(f,gradhess,x1,x2,options); + + //Calculating the values for the 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; + lambda = struct("lower", zl,"upper",zu); + + //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; + lambda = struct("lower",[],"upper",[]); + 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 diff --git a/macros/fmincon.bin b/macros/fmincon.bin Binary files differnew file mode 100644 index 0000000..3ea45b3 --- /dev/null +++ b/macros/fmincon.bin diff --git a/macros/fmincon.sci b/macros/fmincon.sci new file mode 100644 index 0000000..2630683 --- /dev/null +++ b/macros/fmincon.sci @@ -0,0 +1,928 @@ +// 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,lambda,gradient,hessian] = fmincon (varargin) + // Solves a multi-variable constrainted optimization problem + // + // Calling Sequence + // xopt = fmincon(f,x0,A,b) + // xopt = fmincon(f,x0,A,b,Aeq,beq) + // xopt = fmincon(f,x0,A,b,Aeq,beq,lb,ub) + // xopt = fmincon(f,x0,A,b,Aeq,beq,lb,ub,nlc) + // xopt = fmincon(f,x0,A,b,Aeq,beq,lb,ub,nlc,options) + // [xopt,fopt] = fmincon(.....) + // [xopt,fopt,exitflag]= fmincon(.....) + // [xopt,fopt,exitflag,output]= fmincon(.....) + // [xopt,fopt,exitflag,output,lambda]=fmincon(.....) + // [xopt,fopt,exitflag,output,lambda,gradient]=fmincon(.....) + // [xopt,fopt,exitflag,output,lambda,gradient,hessian]=fmincon(.....) + // + // Parameters + // f : a function, representing the objective function of the problem + // x0 : a vector of doubles, containing the starting values of variables of size (1 X n) or (n X 1) where 'n' is the number of Variables + // A : a matrix of doubles, containing the coefficients of linear inequality constraints of size (m X n) where 'm' is the number of linear inequality constraints + // b : a vector of doubles, related to 'A' and containing the the Right hand side equation of the linear inequality constraints of size (m X 1) + // Aeq : a matrix of doubles, containing the coefficients of linear equality constraints of size (m1 X n) where 'm1' is the number of linear equality constraints + // beq : a vector of doubles, related to 'Aeq' and containing the the Right hand side equation of the linear equality constraints of size (m1 X 1) + // lb : a vector of doubles, containing the lower bounds of the variables of size (1 X n) or (n X 1) where 'n' is the number of Variables + // ub : a vector of doubles, containing the upper bounds of the variables of size (1 X n) or (n X 1) where 'n' is the number of Variables + // nlc : a function, representing the Non-linear Constraints functions(both Equality and Inequality) of the problem. It is declared in such a way that non-linear inequality constraints are defined first as a single row vector (c), followed by non-linear equality constraints as another single row vector (ceq). Refer Example for definition of Constraint function. + // options : a list, containing the option for user to specify. See below for details. + // xopt : a vector of doubles, cointating the computed solution of the optimization problem + // fopt : a scalar of double, containing the 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. + // lambda : a structure, containing the Lagrange multipliers of lower bound, upper bound and constraints at the optimized point. See below for details. + // gradient : a vector of doubles, containing the Objective's gradient of the solution. + // hessian : a matrix of doubles, containing the Lagrangian's hessian of the solution. + // + // Description + // Search the minimum of a constrained optimization problem specified by : + // Find the minimum of f(x) such that + // + // <latex> + // \begin{eqnarray} + // &\mbox{min}_{x} + // & f(x) \\ + // & \text{subject to} & A*x \leq b \\ + // & & Aeq*x \ = beq\\ + // & & c(x) \leq 0\\ + // & & ceq(x) \ = 0\\ + // & & lb \leq x \leq ub \\ + // \end{eqnarray} + // </latex> + // + // The routine calls Ipopt for solving the 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", [---], "GradObj", ---, "Hessian", ---, "GradCon", ---);</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>GradObj : a function, representing the gradient function of the Objective in Vector Form.</listitem> + // <listitem>Hessian : a function, representing the hessian function of the Lagrange in Symmetric Matrix Form with Input parameters x, Objective factor and Lambda. Refer Example for definition of Lagrangian Hessian function.</listitem> + // <listitem>GradCon : a function, representing the gradient of the Non-Linear Constraints (both Equality and Inequality) of the problem. It is declared in such a way that gradient of non-linear inequality constraints are defined first as a separate Matrix (cg of size m2 X n or as an empty), followed by gradient of non-linear equality constraints as a separate Matrix (ceqg of size m2 X n or as an empty) where m2 & m3 are number of non-linear inequality and equality constraints respectively.</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 amount of 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> + // + // The lambda data structure contains the Lagrange multipliers at the end + // of optimization. In the current version the values are returned only when the the solution is optimal. + // It has type "struct" and contains the following fields. + // <itemizedlist> + // <listitem>lambda.lower: The Lagrange multipliers for the lower bound constraints.</listitem> + // <listitem>lambda.upper: The Lagrange multipliers for the upper bound constraints.</listitem> + // <listitem>lambda.eqlin: The Lagrange multipliers for the linear equality constraints.</listitem> + // <listitem>lambda.ineqlin: The Lagrange multipliers for the linear inequality constraints.</listitem> + // <listitem>lambda.eqnonlin: The Lagrange multipliers for the non-linear equality constraints.</listitem> + // <listitem>lambda.ineqnonlin: The Lagrange multipliers for the non-linear inequality constraints.</listitem> + // </itemizedlist> + // + // Examples + // //Find x in R^2 such that it minimizes: + // //f(x)= -x1 -x2/3 + // //x0=[0,0] + // //constraint-1 (c1): x1 + x2 <= 2 + // //constraint-2 (c2): x1 + x2/4 <= 1 + // //constraint-3 (c3): x1 - x2 <= 2 + // //constraint-4 (c4): -x1/4 - x2 <= 1 + // //constraint-5 (c5): -x1 - x2 <= -1 + // //constraint-6 (c6): -x1 + x2 <= 2 + // //constraint-7 (c7): x1 + x2 = 2 + // //Objective function to be minimised + // function y=f(x) + // y=-x(1)-x(2)/3; + // endfunction + // //Starting point, linear constraints and variable bounds + // x0=[0 , 0]; + // A=[1,1 ; 1,1/4 ; 1,-1 ; -1/4,-1 ; -1,-1 ; -1,1]; + // b=[2;1;2;1;-1;2]; + // Aeq=[1,1]; + // beq=[2]; + // lb=[]; + // ub=[]; + // nlc=[]; + // //Gradient of objective function + // function y= fGrad(x) + // y= [-1,-1/3]; + // endfunction + // //Hessian of lagrangian + // function y= lHess(x,obj,lambda) + // y= obj*[0,0;0,0] + // endfunction + // //Options + // options=list("GradObj", fGrad, "Hessian", lHess); + // //Calling Ipopt + // [x,fval,exitflag,output,lambda,grad,hessian] =fmincon(f, x0,A,b,Aeq,beq,lb,ub,nlc,options) + // // Press ENTER to continue + // + // Examples + // //Find x in R^3 such that it minimizes: + // //f(x)= x1*x2 + x2*x3 + // //x0=[0.1 , 0.1 , 0.1] + // //constraint-1 (c1): x1^2 - x2^2 + x3^2 <= 2 + // //constraint-2 (c2): x1^2 + x2^2 + x3^2 <= 10 + // //Objective function to be minimised + // function y=f(x) + // y=x(1)*x(2)+x(2)*x(3); + // endfunction + // //Starting point, linear constraints and variable bounds + // x0=[0.1 , 0.1 , 0.1]; + // A=[]; + // b=[]; + // Aeq=[]; + // beq=[]; + // lb=[]; + // ub=[]; + // //Nonlinear constraints + // function [c,ceq]=nlc(x) + // c = [x(1)^2 - x(2)^2 + x(3)^2 - 2 , x(1)^2 + x(2)^2 + x(3)^2 - 10]; + // ceq = []; + // endfunction + // //Gradient of objective function + // function y= fGrad(x) + // y= [x(2),x(1)+x(3),x(2)]; + // endfunction + // //Hessian of the Lagrange Function + // function y= lHess(x,obj,lambda) + // y= obj*[0,1,0;1,0,1;0,1,0] + lambda(1)*[2,0,0;0,-2,0;0,0,2] + lambda(2)*[2,0,0;0,2,0;0,0,2] + // endfunction + // //Gradient of Non-Linear Constraints + // function [cg,ceqg] = cGrad(x) + // cg=[2*x(1) , -2*x(2) , 2*x(3) ; 2*x(1) , 2*x(2) , 2*x(3)]; + // ceqg=[]; + // endfunction + // //Options + // options=list("MaxIter", [1500], "CpuTime", [500], "GradObj", fGrad, "Hessian", lHess,"GradCon", cGrad); + // //Calling Ipopt + // [x,fval,exitflag,output] =fmincon(f, x0,A,b,Aeq,beq,lb,ub,nlc,options) + // // Press ENTER to continue + // + // Examples + // //The below problem is an unbounded problem: + // //Find x in R^3 such that it minimizes: + // //f(x)= -(x1^2 + x2^2 + x3^2) + // //x0=[0.1 , 0.1 , 0.1] + // // x1 <= 0 + // // x2 <= 0 + // // x3 <= 0 + // //Objective function to be minimised + // function y=f(x) + // y=-(x(1)^2+x(2)^2+x(3)^2); + // endfunction + // //Starting point, linear constraints and variable bounds + // x0=[0.1 , 0.1 , 0.1]; + // A=[]; + // b=[]; + // Aeq=[]; + // beq=[]; + // lb=[]; + // ub=[0,0,0]; + // //Options + // options=list("MaxIter", [1500], "CpuTime", [500]); + // //Calling Ipopt + // [x,fval,exitflag,output,lambda,grad,hessian] =fmincon(f, x0,A,b,Aeq,beq,lb,ub,[],options) + // // Press ENTER to continue + // + // Examples + // //The below problem is an infeasible problem: + // //Find x in R^3 such that in minimizes: + // //f(x)=x1*x2 + x2*x3 + // //x0=[1,1,1] + // //constraint-1 (c1): x1^2 <= 1 + // //constraint-2 (c2): x1^2 + x2^2 <= 1 + // //constraint-3 (c3): x3^2 <= 1 + // //constraint-4 (c4): x1^3 = 0.5 + // //constraint-5 (c5): x2^2 + x3^2 = 0.75 + // // 0 <= x1 <=0.6 + // // 0.2 <= x2 <= inf + // // -inf <= x3 <= 1 + // //Objective function to be minimised + // function y=f(x) + // y=x(1)*x(2)+x(2)*x(3); + // endfunction + // //Starting point, linear constraints and variable bounds + // x0=[1,1,1]; + // A=[]; + // b=[]; + // Aeq=[]; + // beq=[]; + // lb=[0 0.2,-%inf]; + // ub=[0.6 %inf,1]; + // //Nonlinear constraints + // function [c,ceq]=nlc(x) + // c=[x(1)^2-1,x(1)^2+x(2)^2-1,x(3)^2-1]; + // ceq=[x(1)^3-0.5,x(2)^2+x(3)^2-0.75]; + // endfunction + // //Gradient of objective function + // function y= fGrad(x) + // y= [x(2),x(1)+x(3),x(2)]; + // endfunction + // //Hessian of the Lagrange Function + // function y= lHess(x,obj,lambda) + // y= obj*[0,1,0;1,0,1;0,1,0] + lambda(1)*[2,0,0;0,0,0;0,0,0] + lambda(2)*[2,0,0;0,2,0;0,0,0] +lambda(3)*[0,0,0;0,0,0;0,0,2] + lambda(4)*[6*x(1 ),0,0;0,0,0;0,0,0] + lambda(5)*[0,0,0;0,2,0;0,0,2]; + // endfunction + // //Gradient of Non-Linear Constraints + // function [cg,ceqg] = cGrad(x) + // cg = [2*x(1),0,0;2*x(1),2*x(2),0;0,0,2*x(3)]; + // ceqg = [3*x(1)^2,0,0;0,2*x(2),2*x(3)]; + // endfunction + // //Options + // options=list("MaxIter", [1500], "CpuTime", [500], "GradObj", fGrad, "Hessian", lHess,"GradCon", cGrad); + // //Calling Ipopt + // [x,fval,exitflag,output,lambda,grad,hessian] =fmincon(f, x0,A,b,Aeq,beq,lb,ub,nlc,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<4 | rhs>13 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while it should be 4,6,8,9,10,11,12,13"), "fmincon", rhs); + error(errmsg) + end + + if (rhs==5 | rhs==7) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while it should be 4,6,8,9,10,11,12,13"), "fmincon", rhs); + error(errmsg) + end + + //Storing the Input Parameters + fun = varargin(1); + x0 = varargin(2); + A = varargin(3); + b = varargin(4); + Aeq = []; + beq = []; + lb = []; + ub = []; + nlc = []; + + if (rhs>4) then + Aeq = varargin(5); + beq = varargin(6); + end + + if (rhs>6) then + lb = varargin(7); + ub = varargin(8); + end + + if (rhs>8) then + nlc = varargin(9); + end + + //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 (1st Parameter)"), "fmincon"); + 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 (2nd Parameter)"), "fmincon"); + 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 (Starting Point) or Starting Point cannot be Empty"), "fmincon"); + error(errmsg); + end + + 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); + + //To check the match between fun (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"), "fmincon"); + 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 the 3rd Input argument (A) is a Matrix/Vector + if (type(A) ~= 1) then + errmsg = msprintf(gettext("%s: Expected Matrix/Vector for Constraint Matrix A (3rd parameter)"), "fmincon"); + error(errmsg); + end + + //To check for correct size of A(3rd paramter) + if(size(A,2)~=s(2) & size(A,2)~=0) then + errmsg = msprintf(gettext("%s: Expected Matrix of size (No of linear inequality constraints X No of Variables) or an Empty Matrix for Linear Inequality Constraint coefficient Matrix A"), "fmincon"); + error(errmsg); + end + + s1=size(A); + + //To check whether the 4th Input argument (b) is a vector/scalar + if (type(b) ~= 1) then + errmsg = msprintf(gettext("%s: Expected Vector/Scalar for b (4th Parameter)"), "fmincon"); + error(errmsg); + end + + //To check for the correct size of b (4th paramter) and convert it into a column vector which is required for Ipopt + if(s1(2)==0) then + if(size(b,2)~=0) then + errmsg = msprintf(gettext("%s: As Linear Inequality Constraint coefficient Matrix A (3rd parameter) is empty, b (4th Parameter) should also be empty"), "fmincon"); + error(errmsg); + end + else + if((size(b,1)~=1) & (size(b,2)~=1)) then + errmsg = msprintf(gettext("%s: Expected Non empty Row/Column Vector for b (4th Parameter) for your Inputs "), "fmincon"); + error(errmsg); + elseif(size(b,1)~=s1(1) & size(b,2)==1) then + errmsg = msprintf(gettext("%s: Expected Column Vector (number of linear inequality constraints X 1) for b (4th Parameter) "), "fmincon"); + error(errmsg); + elseif(size(b,1)==s1(1) & size(b,2)==1) then + b=b; + elseif(size(b,1)==1 & size(b,2)~=s1(1)) then + errmsg = msprintf(gettext("%s: Expected Row Vector (1 X number of linear inequality constraints) for b (4th Parameter) "), "fmincon"); + error(errmsg); + elseif(size(b,1)==1 & size(b,2)==s1(1)) then + b=b'; + end + end + + //To check whether the 5th Input argument (Aeq) is a matrix/vector + if (type(Aeq) ~= 1) then + errmsg = msprintf(gettext("%s: Expected Matrix/Vector for Equality Constraint Matrix Aeq (5th Parameter)"), "fmincon"); + error(errmsg); + end + + //To check for the correct size of Aeq (5th paramter) + if(size(Aeq,2)~=s(2) & size(Aeq,2)~=0) then + errmsg = msprintf(gettext("%s: Expected Matrix of size (No of linear equality constraints X No of Variables) or an Empty Matrix for Linear Equality Constraint coefficient Matrix Aeq"), "fmincon"); + error(errmsg); + end + + s2=size(Aeq); + + //To check whether the 6th Input argument(beq) is a vector/scalar + if (type(beq) ~= 1) then + errmsg = msprintf(gettext("%s: Expected Vector/Scalar for beq (6th Parameter)"), "fmincon"); + error(errmsg); + end + + //To check for the correct size of beq(6th paramter) and convert it into a column vector which is required for Ipopt + if(s2(2)==0) then + if(size(beq,2)~=0) then + errmsg = msprintf(gettext("%s: As Linear Equality Constraint coefficient Matrix Aeq (5th parameter) is empty, beq (6th Parameter) should also be empty"), "fmincon"); + error(errmsg); + end + else + if((size(beq,1)~=1) & (size(beq,2)~=1)) then + errmsg = msprintf(gettext("%s: Expected Non empty Row/Column Vector for beq (6th Parameter)"), "fmincon"); + error(errmsg); + elseif(size(beq,1)~=s2(1) & size(beq,2)==1) then + errmsg = msprintf(gettext("%s: Expected Column Vector (number of linear equality constraints X 1) for beq (6th Parameter) "), "fmincon"); + error(errmsg); + elseif(size(beq,1)==s2(1) & size(beq,2)==1) then + beq=beq; + elseif(size(beq,1)==1 & size(beq,2)~=s2(1)) then + errmsg = msprintf(gettext("%s: Expected Row Vector (1 X number of linear equality constraints) for beq (6th Parameter) "), "fmincon"); + error(errmsg); + elseif(size(beq,1)==1 & size(beq,2)==s2(1)) then + beq=beq'; + end + end + + + //To check whether the 7th Input argument (lb) is a vector/scalar + if (type(lb) ~= 1) then + errmsg = msprintf(gettext("%s: Expected Vector/Scalar for Lower Bound Vector (7th Parameter)"), "fmincon"); + error(errmsg); + end + + //To check for the correct size and data of lb (7th paramter) and convert it into a column vector as required by Ipopt + if (size(lb,2)==0) then + lb = repmat(-%inf,1,s(2)); + end + + if (size(lb,1)~=1) & (size(lb,2)~=1) then + errmsg = msprintf(gettext("%s: Lower Bound (7th Parameter) should be a vector"), "fmincon"); + error(errmsg); + elseif(size(lb,1)~=s(2) & size(lb,2)==1) then + errmsg = msprintf(gettext("%s: Expected Column Vector (number of Variables X 1) for lower bound (7th Parameter) "), "fmincon"); + error(errmsg); + elseif(size(lb,1)==s(2) & size(lb,2)==1) then + lb=lb; + elseif(size(lb,1)==1 & size(lb,2)~=s(2)) then + errmsg = msprintf(gettext("%s: Expected Row Vector (1 X number of Variables) for lower bound (7th Parameter) "), "fmincon"); + error(errmsg); + elseif(size(lb,1)==1 & size(lb,2)==s(2)) then + lb=lb'; + end + + //To check whether the 8th Input argument (ub) is a vector/scalar + if (type(ub) ~= 1) then + errmsg = msprintf(gettext("%s: Expected Vector/Scalar for Upper Bound Vector (8th Parameter)"), "fmincon"); + error(errmsg); + end + + //To check for the correct size and data of ub (8th paramter) and convert it into a column vector as required by Ipopt + if (size(ub,2)==0) then + ub = repmat(%inf,1,s(2)); + end + + if (size(ub,1)~=1)& (size(ub,2)~=1) then + errmsg = msprintf(gettext("%s: Upper Bound (8th Parameter) should be a vector"), "fmincon"); + error(errmsg); + elseif(size(ub,1)~=s(2) & size(ub,2)==1) then + errmsg = msprintf(gettext("%s: Expected Column Vector (number of Variables X 1) for upper bound (8th Parameter) "), "fmincon"); + error(errmsg); + elseif(size(ub,1)==s(2) & size(ub,2)==1) then + ub=ub; + elseif(size(ub,1)==1 & size(ub,2)~=s(2)) then + errmsg = msprintf(gettext("%s: Expected Row Vector (1 X number of Variables) for upper bound (8th Parameter) "), "fmincon"); + error(errmsg); + elseif(size(ub,1)==1 & size(ub,2)==s(2)) then + ub=ub'; + end + + //To check the contents of lb & ub (7th & 8th Parameter) + for i = 1:s(2) + if (lb(i) == %inf) then + errmsg = msprintf(gettext("%s: Value of Lower Bound can not be infinity"), "fmincon"); + error(errmsg); + end + + if (ub(i) == -%inf) then + errmsg = msprintf(gettext("%s: Value of Upper Bound can not be negative infinity"), "fmincon"); + error(errmsg); + end + + if(ub(i)-lb(i)<=1e-6) then + errmsg = msprintf(gettext("%s: Difference between Upper Bound and Lower bound should be atleast > 10^6 for variable number= %d "), "fmincon", i); + error(errmsg) + end + end + + //To check whether the 10th Input argument (nlc) is a function or an empty matrix + if (type(nlc) == 1 & size(nlc,2)==0 ) then + addnlc=[]; + addnlc1=[]; + no_nlc=0; + no_nlic=0; + no_nlec=0; + elseif (type(nlc) == 13 | type(nlc) == 11) then + + if(execstr('[sample_c,sample_ceq] = nlc(x0)','errcatch')==21) + errmsg = msprintf(gettext("%s: Non-Linear Constraint function(9th Parameter) and x0(2nd Parameter) did not match"), "fmincon"); + error(errmsg); + end + [sample_c,sample_ceq] = nlc(x0); + + if (size(sample_c,1)~=1 & size(sample_c,1)~=0) then + errmsg = msprintf(gettext("%s: Definition of c in Non-Linear Constraint function(9th Parameter) should be in the form of Row Vector or Empty Vector"), "fmincon"); + error(errmsg) + end + + if (size(sample_ceq,1)~=1 & size(sample_ceq,1)~=0) then + errmsg = msprintf(gettext("%s: Definition of ceq in Non-Linear Constraint function(9th Parameter) should be in the form of Row Vector or Empty Vector"), "fmincon"); + error(errmsg) + end + + no_nlic = size(sample_c,2); + no_nlec = size(sample_ceq,2); + no_nlc = no_nlic + no_nlec; + + //Constructing a single output variable function for nlc + function y = addnlc(x) + [c,ceq] = nlc(x); + y = [c';ceq']; + endfunction + + //To check the addnlc function + if(execstr('sample_allcon = addnlc(x0)','errcatch')==21) + errmsg = msprintf(gettext("%s: Non-Linear Constraint function(9th Parameter) and x0(2nd Parameter) did not match"), "fmincon"); + error(errmsg); + end + sample_allcon = addnlc(x0); + + if (size(sample_allcon,1)==0 & size(sample_allcon,2)==0) then + + elseif (size(sample_allcon,1)~=no_nlc | size(sample_allcon,2)~=1) then + errmsg = msprintf(gettext("%s: Please check the Non-Linear Constraint function (9th Parameter) function"), "fmincon"); + error(errmsg) + end + + //Constructing a nlc function with error deduction + function [y,check] = addnlc1(x) + if(execstr('y = addnlc(x)','errcatch')==32 | execstr('y = addnlc(x)','errcatch')==27) + y = 0; + check=1; + else + y = addnlc(x); + if (isreal(y)==%F) then + y = 0; + check=1; + else + check=0; + end + end + endfunction + + else + errmsg = msprintf(gettext("%s: Non Linear Constraint (9th Parameter) should be a function or an Empty Matrix"), "fmincon"); + error(errmsg) + end + + //To check whether options has been entered by the user + if ( rhs<10 ) then + param = list(); + else + param =varargin(10); //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: Options (10th parameter) should be a list"), "fmincon"); + 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 Options (list) should be even"), "fmincon"); + 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 + elseif t==2 then //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 + elseif t==3 then //To return Gradient + if(execstr('y=numderivative(addnlc,x)','errcatch')==10000) + y=0; + check=1; + else + y=numderivative(addnlc,x); + if (isreal(y)==%F) then + y=0; + check=1; + else + check=0; + end + end + elseif t==4 then //To return Hessian + if(execstr('[grad,y]=numderivative(addnlc,x)','errcatch')==10000) + y=0; + check=1; + else + [grad,y]=numderivative(addnlc,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; + flag3=0; + + //Function for Gradient and Hessian + fGrad=[]; + fGrad1=[]; + lHess=[]; + lHess1=[]; + cGrad=[]; + addcGrad=[]; + addcGrad1=[]; + + //To check the user entry for options and storing 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 "GradObj" then + flag1=1; + fGrad=param(2*i); + case "Hessian" then + flag2=1; + lHess=param(2*i); + case "GradCon" then + flag3=1; + cGrad=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) ~= 11 & type(fGrad) ~= 13) then + errmsg = msprintf(gettext("%s: Expected function for Gradient of Objective"), "fmincon"); + error(errmsg); + end + + if(execstr('sample_fGrad=fGrad(x0)','errcatch')==21) + errmsg = msprintf(gettext("%s: Gradient function of Objective and x0 did not match "), "fmincon"); + error(errmsg); + end + + sample_fGrad=fGrad(x0); + + if (size(sample_fGrad,1)==s(2) & size(sample_fGrad,2)==1) then + elseif (size(sample_fGrad,1)==1 & size(sample_fGrad,2)==s(2)) then + elseif (size(sample_fGrad,1)~=1 & size(sample_fGrad,2)~=1) then + errmsg = msprintf(gettext("%s: Wrong Input for Objective Gradient function(10th Parameter)---->Vector function is Expected"), "fmincon"); + 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(lHess) ~= 11 & type(lHess) ~= 13) then + errmsg = msprintf(gettext("%s: Expected function for Hessian of Objective"), "fmincon"); + error(errmsg); + end + if(execstr('sample_lHess=lHess(x0,1,1:no_nlc)','errcatch')==21) + errmsg = msprintf(gettext("%s: Hessian function of Objective and x0 did not match "), "fmincon"); + error(errmsg); + end + sample_lHess=lHess(x0,1,1:no_nlc); + if(size(sample_lHess,1)~=s(2) | size(sample_lHess,2)~=s(2)) then + errmsg = msprintf(gettext("%s: Wrong Input for Objective Hessian function(10th Parameter)---->Symmetric Matrix function is Expected "), "fmincon"); + error(errmsg); + end + + function [y,check] = lHess1(x,obj,lambda) + if(execstr('y=lHess(x,obj,lambda)','errcatch')==32 | execstr('y=lHess(x,obj,lambda)','errcatch')==27) + y = 0; + check=1; + else + y=lHess(x,obj,lambda); + if (isreal(y)==%F) then + y = 0; + check=1; + else + check=0; + end + end + endfunction + + end + if (flag3==1) then + if (type(cGrad) ~= 11 & type(cGrad) ~= 13) then + errmsg = msprintf(gettext("%s: Expected function for Gradient of Constraint function"), "fmincon"); + error(errmsg); + end + + if(execstr('[sample_cGrad,sample_ceqg]=cGrad(x0)','errcatch')==21) + errmsg = msprintf(gettext("%s: Gradient function of Constraint and x0 did not match "), "fmincon"); + error(errmsg); + end + [sample_cGrad,sample_ceqg]=cGrad(x0); + + if (size(sample_cGrad,2)==0) then + elseif (size(sample_cGrad,1)~=no_nlic | size(sample_cGrad,2)~=s(2)) then + errmsg = msprintf(gettext("%s: Definition of (cGrad) in Non-Linear Constraint function(10th Parameter) should be in the form of (m X n) or Empty Matrix where m is number of Non- linear inequality constraints and n is number of Variables"), "fmincon"); + error(errmsg); + end + + if (size(sample_ceqg,2)==0) then + elseif (size(sample_ceqg,1)~=no_nlec | size(sample_ceqg,2)~=s(2)) then + errmsg = msprintf(gettext("%s: Definition of (ceqg) in Non-Linear Constraint function(10th Parameter) should be in the form of (m X n) or Empty Matrix where m is number of Non- linear equality constraints and n is number of Variables"), "fmincon"); + error(errmsg); + end + + function y = addcGrad(x) + [sample_cGrad,sample_ceqg] = cGrad(x); + y = [sample_cGrad;sample_ceqg]; + endfunction + + sample_addcGrad=addcGrad(x0); + if(size(sample_addcGrad,1)~=no_nlc | size(sample_addcGrad,2)~=s(2)) then + errmsg = msprintf(gettext("%s: Wrong Input for Constraint Gradient function(10th Parameter) (Refer Help)"), "fmincon"); + error(errmsg); + end + + function [y,check] = addcGrad1(x) + if(execstr('y=addcGrad(x)','errcatch')==32 | execstr('y=addcGrad(x)','errcatch')==27) + y = 0; + check=1; + else + y=addcGrad(x); + if (isreal(y)==%F) then + y = 0; + check=1; + else + check=0; + end + end + endfunction + end + + //To Convert the Gradient and Hessian into Error Debugable form + + + //Dummy variable which is used by Ipopt + empty=0; + + //Calling the Ipopt function for solving the above problem + [xopt,fopt,status,iter,cpu,obj_eval,dual,lambda1,zl,zu,gradient,hessian1] = solveminconp (f,gradhess,A,b,Aeq,beq,lb,ub,no_nlc,no_nlic,addnlc1,flag1,fGrad1,flag2,lHess1,flag3,addcGrad1,x0,options,empty) + + //Calculating the values for the 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; + lambda = struct("lower", zl,"upper",zu,"ineqlin",[],"eqlin",[],"ineqnonlin",[],"eqnonlin",[]); + + if (no_nlic ~= 0) then + for i = 1:no_nlic + lambda.ineqnonlin (i) = lambda1(i) + end + lambda.ineqnonlin = lambda.ineqnonlin' + end + + if (no_nlec ~= 0) then + j=1; + for i = no_nlic+1 : no_nlc + lambda.eqnonlin (j) = lambda1(i) + j= j+1; + end + lambda.eqnonlin = lambda.eqnonlin' + end + + if (Aeq ~=[]) then + j=1; + for i = no_nlc+1 : no_nlc + size(Aeq,1) + lambda.eqlin (j) = lambda1(i) + j= j+1; + end + lambda.eqlin = lambda.eqlin' + end + + if (A ~=[]) then + j=1; + for i = no_nlc+ size(Aeq,1)+ 1 : no_nlc + size(Aeq,1) + size(A,1) + lambda.ineqlin (j) = lambda1(i) + j= j+1; + end + lambda.ineqlin = lambda.ineqlin' + end + + //Converting hessian of order (1 x (numberOfVariables)^2) received from Ipopt to order (numberOfVariables x numberOfVariables) + s1=size(gradient) + for i =1:s1(2) + for j =1:s1(2) + hessian(i,j)= hessian1(j+((i-1)*s1(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; + lambda = struct("lower",[],"upper",[],"ineqlin",[],"eqlin",[],"ineqnonlin",[],"eqnonlin",[]); + 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 11 then + printf("\nProblem has too few degrees of freedom.\n"); + case 12 then + printf("\nInvalid option thrown back by Ipopt\n"); + case 13 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 diff --git a/macros/fminimax.bin b/macros/fminimax.bin Binary files differnew file mode 100644 index 0000000..c40fcce --- /dev/null +++ b/macros/fminimax.bin diff --git a/macros/fminimax.sci b/macros/fminimax.sci new file mode 100644 index 0000000..5084475 --- /dev/null +++ b/macros/fminimax.sci @@ -0,0 +1,511 @@ +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2015 - IIT Bombay - FOSSEE +// +// Authors: Animesh Baranawal +// Organization: FOSSEE, IIT Bombay +// Email: animeshbaranawal@gmail.com +// 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 + +function [x,fval,maxfval,exitflag,output,lambda] = fminimax(varargin) + // Solves minimax constraint problem + // + // Calling Sequence + // x = fminimax(fun,x0) + // x = fminimax(fun,x0,A,b) + // x = fminimax(fun,x0,A,b,Aeq,beq) + // x = fminimax(fun,x0,A,b,Aeq,beq,lb,ub) + // x = fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlinfun) + // x = fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlinfun,options) + // [x, fval] = fmincon(.....) + // [x, fval, maxfval]= fmincon(.....) + // [x, fval, maxfval, exitflag]= fmincon(.....) + // [x, fval, maxfval, exitflag, output]= fmincon(.....) + // [x, fval, maxfval, exitflag, output, lambda]= fmincon(.....) + // + // Parameters + // fun: The function to be minimized. fun is a function that accepts a vector x and returns a vector F, the objective functions evaluated at x. + // x0: a nx1 or 1xn matrix of doubles, where n is the number of variables, the initial guess for the optimization algorithm + // A: a nil x n matrix of doubles, where n is the number of variables and nil is the number of linear inequalities. If A==[] and b==[], it is assumed that there is no linear inequality constraints. If (A==[] & b<>[]), fminimax generates an error (the same happens if (A<>[] & b==[])) + // b: a nil x 1 matrix of doubles, where nil is the number of linear inequalities + // Aeq: a nel x n matrix of doubles, where n is the number of variables and nel is the number of linear equalities. If Aeq==[] and beq==[], it is assumed that there is no linear equality constraints. If (Aeq==[] & beq<>[]), fminimax generates an error (the same happens if (Aeq<>[] & beq==[])) + // beq: a nel x 1 matrix of doubles, where nel is the number of linear equalities + // lb: a nx1 or 1xn matrix of doubles, where n is the number of variables. The lower bound for x. If lb==[], then the lower bound is automatically set to -inf + // ub: a nx1 or 1xn matrix of doubles, where n is the number of variables. The upper bound for x. If ub==[], then the upper bound is automatically set to +inf + // nonlinfun: function that computes the nonlinear inequality constraints c(x) <= 0 and nonlinear equality constraints ceq(x) = 0. + // x: a nx1 matrix of doubles, the computed solution of the optimization problem + // fval: a vector of doubles, the value of fun at x + // maxfval: a 1x1 matrix of doubles, the maximum value in vector fval + // exitflag: a 1x1 matrix of floating point integers, the exit status + // output: a struct, the details of the optimization process + // lambda: a struct, the Lagrange multipliers at optimum + // options: a list, containing the option for user to specify. See below for details. + // + // Description + // fminimax minimizes the worst-case (largest) value of a set of multivariable functions, starting at an initial estimate. This is generally referred to as the minimax problem. + // + // <latex> + // \min_{x} \max_{i} F_{i}(x)\: \textrm{such that} \:\begin{cases} + // & c(x) \leq 0 \\ + // & ceq(x) = 0 \\ + // & A.x \leq b \\ + // & Aeq.x = beq \\ + // & minmaxLb \leq x \leq minmaxUb + // \end{cases} + // </latex> + // + // Currently, fminimax calls fmincon which uses the ip-opt algorithm. + // + // max-min problems can also be solved with fminimax, using the identity + // + // <latex> + // \max_{x} \min_{i} F_{i}(x) = -\min_{x} \max_{i} \left( -F_{i}(x) \right) + // </latex> + // + // 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", [---], "GradObj", ---, "GradCon", ---);</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>GradObj : a function, representing the gradient function of the Objective in Vector Form.</listitem> + // <listitem>GradCon : a function, representing the gradient of the Non-Linear Constraints (both Equality and Inequality) of the problem. It is declared in such a way that gradient of non-linear inequality constraints are defined first as a separate Matrix (cg of size m2 X n or as an empty), followed by gradient of non-linear equality constraints as a separate Matrix (ceqg of size m2 X n or as an empty) where m2 & m3 are number of non-linear inequality and equality constraints respectively.</listitem> + // <listitem>Default Values : options = list("MaxIter", [3000], "CpuTime", [600]);</listitem> + // </itemizedlist> + // + // The objective function must have header : + // <programlisting> + // F = fun(x) + // </programlisting> + // where x is a n x 1 matrix of doubles and F is a m x 1 matrix of doubles where m is the total number of objective functions inside F. + // On input, the variable x contains the current point and, on output, the variable F must contain the objective function values. + // + // By default, the gradient options for fminimax are turned off and and fmincon does the gradient opproximation of minmaxObjfun. In case the GradObj option is off and GradConstr option is on, fminimax approximates minmaxObjfun gradient using numderivative toolbox. + // + // If we can provide exact gradients, we should do so since it improves the convergence speed of the optimization algorithm. + // + // Furthermore, we must enable the "GradObj" option with the statement : + // <programlisting> + // minimaxOptions = list("GradObj",fGrad); + // </programlisting> + // This will let fminimax know that the exact gradient of the objective function is known, so that it can change the calling sequence to the objective function. Note that, fGrad should be mentioned in the form of N x n where n is the number of variables, N is the number of functions in objective function. + // + // The constraint function must have header : + // <programlisting> + // [c, ceq] = confun(x) + // </programlisting> + // where x is a n x 1 matrix of dominmaxUbles, c is a 1 x nni matrix of doubles and ceq is a 1 x nne matrix of doubles (nni : number of nonlinear inequality constraints, nne : number of nonlinear equality constraints). + // On input, the variable x contains the current point and, on output, the variable c must contain the nonlinear inequality constraints and ceq must contain the nonlinear equality constraints. + // + // By default, the gradient options for fminimax are turned off and and fmincon does the gradient opproximation of confun. In case the GradObj option is on and GradCons option is off, fminimax approximates confun gradient using numderivative toolbox. + // + // If we can provide exact gradients, we should do so since it improves the convergence speed of the optimization algorithm. + // + // Furthermore, we must enable the "GradCon" option with the statement : + // <programlisting> + // minimaxOptions = list("GradCon",confunGrad); + // </programlisting> + // This will let fminimax know that the exact gradient of the objective function is known, so that it can change the calling sequence to the objective function. + // + // The constraint derivative function must have header : + // <programlisting> + // [dc,dceq] = confungrad(x) + // </programlisting> + // where dc is a nni x n matrix of doubles and dceq is a nne x n matrix of doubles. + // + // 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 amount of 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> + // + // The lambda data structure contains the Lagrange multipliers at the end + // of optimization. In the current version the values are returned only when the the solution is optimal. + // It has type "struct" and contains the following fields. + // <itemizedlist> + // <listitem>lambda.lower: The Lagrange multipliers for the lower bound constraints.</listitem> + // <listitem>lambda.upper: The Lagrange multipliers for the upper bound constraints.</listitem> + // <listitem>lambda.eqlin: The Lagrange multipliers for the linear equality constraints.</listitem> + // <listitem>lambda.ineqlin: The Lagrange multipliers for the linear inequality constraints.</listitem> + // <listitem>lambda.eqnonlin: The Lagrange multipliers for the non-linear equality constraints.</listitem> + // <listitem>lambda.ineqnonlin: The Lagrange multipliers for the non-linear inequality constraints.</listitem> + // </itemizedlist> + // + // Examples + // // A basic case : + // // we provide only the objective function and the nonlinear constraint + // // function + // function f = myfun(x) + // f(1)= 2*x(1)^2 + x(2)^2 - 48*x(1) - 40*x(2) + 304; //Objectives + // f(2)= -x(1)^2 - 3*x(2)^2; + // f(3)= x(1) + 3*x(2) -18; + // f(4)= -x(1) - x(2); + // f(5)= x(1) + x(2) - 8; + // endfunction + // // The initial guess + // x0 = [0.1,0.1]; + // // The expected solution : only 4 digits are guaranteed + // //xopt = [4 4] + // //fopt = [0 -64 -2 -8 0] + // maxfopt = 0 + // // Run fminimax + // [xopt,fopt,maxfval,exitflag,output,lambda] = fminimax(myfun, x0) + // // Press ENTER to continue + // + // Examples + // // A case where we provide the gradient of the objective + // // functions and the Jacobian matrix of the constraints. + // // The objective function and its gradient + // function f = myfun(x) + // f(1)= 2*x(1)^2 + x(2)^2 - 48*x(1) - 40*x(2) + 304; + // f(2)= -x(1)^2 - 3*x(2)^2; + // f(3)= x(1) + 3*x(2) -18; + // f(4)= -x(1) - x(2); + // f(5)= x(1) + x(2) - 8; + // endfunction + // // Defining gradient of myfun + // function G = myfungrad(x) + // G = [ 4*x(1) - 48, -2*x(1), 1, -1, 1; + // 2*x(2) - 40, -6*x(2), 3, -1, 1; ]' + // endfunction + // // The nonlinear constraints and the Jacobian + // // matrix of the constraints + // function [c,ceq] = confun(x) + // // Inequality constraints + // c = [1.5 + x(1)*x(2) - x(1) - x(2), -x(1)*x(2) - 10] + // // No nonlinear equality constraints + // ceq=[] + // endfunction + // // Defining gradient of confungrad + // function [DC,DCeq] = cgrad(x) + // // DC(:,i) = gradient of the i-th constraint + // // DC = [ + // // Dc1/Dx1 Dc1/Dx2 + // // Dc2/Dx1 Dc2/Dx2 + // // ] + // DC= [ + // x(2)-1, -x(2) + // x(1)-1, -x(1) + // ]' + // DCeq = []' + // endfunction + // // Test with both gradient of objective and gradient of constraints + // minimaxOptions = list("GradObj",myfungrad,"GradCon",cgrad); + // // The initial guess + // x0 = [0,10]; + // // The expected solution : only 4 digits are guaranteed + // //xopt = [0.92791 7.93551] + // //fopt = [6.73443 -189.778 6.73443 -8.86342 0.86342] + // maxfopt = 6.73443 + // // Run fminimax + // [xopt,fopt,maxfval,exitflag,output] = fminimax(myfun,x0,[],[],[],[],[],[], confun, minimaxOptions) + // Authors + // Animesh Baranawal + + // Check number of input and output arguments + [minmaxLhs,minmaxRhs] = argn() + fminimaxCheckrhs("fminimax", minmaxRhs, [2 4 6 8 9 10]) + fminimaxChecklhs("fminimax", minmaxLhs, 1:7) + + // Proper initialisation of objective function + minmaxObjfun = varargin(1) + fminimaxChecktype("fminimax", minmaxObjfun, "minmaxObjfun", 1, "function") + + // Proper initialisation of starting point + minmaxStartpoint = varargin(2) + fminimaxChecktype("fminimax", minmaxStartpoint, "minmaxStartpoint", 2, "constant") + + minmaxNumvar = size(minmaxStartpoint,"*") + fminimaxCheckvector("fminimax", minmaxStartpoint, "minmaxStartpoint", 2, minmaxNumvar) + minmaxStartpoint = minmaxStartpoint(:) + + // Proper initialisation of A and b + if(minmaxRhs < 3) then // if A and b are not provided, declare as empty + minmaxA = [] + minmaxB = [] + else + minmaxA = varargin(3) + minmaxB = varargin(4) + end + + fminimaxChecktype("fminimax", minmaxA, "A", 3, "constant") + fminimaxChecktype("fminimax", minmaxB, "b", 4, "constant") + + // Check if A and b of proper dimensions + if(minmaxA <> [] & minmaxB == []) then + errmsg = msprintf(gettext("%s: Incompatible input arguments #%d and #%d: matrix A is empty, but the column vector b is not empty"), "fminimax", 3, 4) + error(errmsg) + end + + if(minmaxA == [] & minmaxB <> []) then + errmsg = msprintf(gettext("%s: Incompatible input arguments #%d and #%d: matrix A is not empty, but the column vector b is empty"), "fminimax", 3, 4) + error(errmsg) + end + + minmaxNumrowA = size(minmaxA,"r") + if(minmaxA <> []) then + fminimaxCheckdims("fminimax", minmaxA, "A", 3, [minmaxNumrowA minmaxNumvar]) + fminimaxCheckvector("fminimax", minmaxB, "b", 4, minmaxNumrowA) + minmaxB = minmaxB(:) + end + + // Proper initialisation of Aeq and beq + if(minmaxRhs < 5) then // if Aeq and beq are not provided, declare as empty + minmaxAeq = [] + minmaxBeq = [] + else + minmaxAeq = varargin(5) + minmaxBeq = varargin(6) + end + + fminimaxChecktype("fminimax", minmaxAeq, "Aeq", 5, "constant") + fminimaxChecktype("fminimax", minmaxBeq, "beq", 6, "constant") + + // Check if Aeq and beq of proper dimensions + if(minmaxAeq <> [] & minmaxBeq == []) then + errmsg = msprintf(gettext("%s: Incompatible input arguments #%d and #%d: matrix Aeq is empty, but the column vector beq is not empty"), "fminimax", 5, 6) + error(errmsg) + end + + if(minmaxAeq == [] & minmaxBeq <> []) then + errmsg = msprintf(gettext("%s: Incompatible input arguments #%d and #%d: matrix Aeq is not empty, but the column vector beq is empty"), "fminimax", 5, 6) + error(errmsg) + end + + minmaxNumrowAeq = size(minmaxAeq,"r") + if(minmaxAeq <> []) then + fminimaxCheckdims("fminimax", minmaxAeq, "Aeq", 5, [minmaxNumrowAeq minmaxNumvar]) + fminimaxCheckvector("fminimax", minmaxBeq, "beq", 6, minmaxNumrowAeq) + minmaxBeq = minmaxBeq(:) + end + + // Proper initialisation of minmaxLb and minmaxUb + if(minmaxRhs < 7) then // if minmaxLb and minmaxUb are not provided, declare as empty + minmaxLb = [] + minmaxUb = [] + else + minmaxLb = varargin(7) + minmaxUb = varargin(8) + end + + fminimaxChecktype("fminimax", minmaxLb, "lb", 7, "constant") + fminimaxChecktype("fminimax", minmaxUb, "ub", 8, "constant") + + // Check dimensions of minmaxLb and minmaxUb + if(minmaxLb <> []) then + fminimaxCheckvector("fminimax", minmaxLb, "lb", 7, minmaxNumvar) + minmaxLb = minmaxLb(:) + end + + if(minmaxUb <> []) then + fminimaxCheckvector("fminimax", minmaxUb, "ub", 8, minmaxNumvar) + minmaxUb = minmaxUb(:) + end + + // Proper Initialisation of minmaxNonlinfun + if(minmaxRhs < 9) then // if minmaxNonlinfun is not provided, declare as empty + minmaxNonlinfun = [] + else + minmaxNonlinfun = varargin(9) + end + + // fmincon library of scilab gives error when 'c' component of minmaxNonlinfun empty + // add a trivial case of -5 <= 0 to c to bypass this error + if(minmaxNonlinfun == []) then + function [c,ceq] = t(z) + c = [] + ceq = [] + endfunction + minmaxNonlinfun = t + end + + fminimaxChecktype("fminimax", minmaxNonlinfun, "nonlinfun", 9, "function") + + //To check, Whether minimaxOptions is been entered by user + if ( minmaxRhs<10 ) then + minmaxUserOptions = list(); + else + minmaxUserOptions = varargin(10); //Storing the 3rd Input minmaxUserOptionseter in intermediate list named 'minmaxUserOptions' + end + + //If minimaxOptions is entered then checking its type for 'list' + if (type(minmaxUserOptions) ~= 15) then + errmsg = msprintf(gettext("%s: minimaxOptions (10th parameter) should be a list"), "fminimax"); + error(errmsg); + end + + //If minimaxOptions is entered then checking whether even number of entires are entered + if (modulo(size(minmaxUserOptions),2)) then + errmsg = msprintf(gettext("%s: Size of minimaxOptions (list) should be even"), "fminimax"); + error(errmsg); + end + + //Flags to check whether Gradient is "ON"/"OFF" and store values of user options + flag1=0; + flag2=0; + minmaxMaxIter = 3000 + minmaxCPU = 600 + minmaxFGrad=[]; + minmaxCGrad=[]; + + //To check the User Entry for Options and storing it + for i = 1:(size(minmaxUserOptions))/2 + select minmaxUserOptions(2*i-1) + case "MaxIter" then + minmaxIter = minmaxUserOptions(2*i); //Setting the Maximum Iteration as per user entry + + case "CpuTime" then + minmaxCPU = minmaxUserOptions(2*i); //Setting the Maximum CPU Time as per user entry + + case "GradObj" then + flag1=1; + minmaxFGrad = minmaxUserOptions(2*i); + + case "GradCon" then + flag2=1; + minmaxCGrad = minmaxUserOptions(2*i); + + else + errmsg = msprintf(gettext("%s: Unrecognized minmaxUserOptionseter name ''%s''."), "fminimax", minmaxUserOptions(2*i-1)); + error(errmsg) + end + end + + // Checking if minmaxFGrad and minmaxCGrad are functions + if (flag1==1) then + if (type(minmaxFGrad) ~= 11 & type(minmaxFGrad) ~= 13) then + errmsg = msprintf(gettext("%s: Expected function for Gradient of Objective"), "fminimax"); + error(errmsg); + end + end + + if (flag2==1) then + if (type(minmaxCGrad) ~= 11 & type(minmaxCGrad) ~= 13) then + errmsg = msprintf(gettext("%s: Expected function for Gradient of Nonlinfun"), "fminimax"); + error(errmsg); + end + end + + + // Reformulating the problem fminimax to fmincon + minmaxObjfunval = minmaxObjfun(minmaxStartpoint) + minmaxStartpoint(minmaxNumvar+1) = max(minmaxObjfunval) + + if(minmaxA <> []) then + minmaxA = [minmaxA, zeros(minmaxNumrowA,1)] + end + if(minmaxAeq <> []) then + minmaxAeq = [minmaxAeq, zeros(minmaxNumrowAeq,1)] + end + if(minmaxLb <> []) then + minmaxLb(minmaxNumvar+1) = -%inf + end + if(minmaxUb <> []) then + minmaxUb(minmaxNumvar+1) = +%inf + end + + // function handle defining the additional inequalities + function temp = minmaxAddIneq(z) + temp = minmaxObjfun(z) - z(minmaxNumvar+1) + endfunction + + // function handle defining new objective function + function newfunc = newObjfun(z) + newfunc = z(minmaxNumvar+1) + endfunction + + // function handle defining add_ineq derivative using numderivative + function func = minmaxObjDer(z) + func = numderivative(minmaxAddIneq,z) + endfunction + + // function handle defining minmaxNonlinfun derivative using numderivative + function [dc,dceq] = minmaxNonlinDer(z) + // function handle extracting c and ceq components from minmaxNonlinfun + function foo = minmaxC(z) + [foo,tmp1] = minmaxNonlinfun(z) + foo = foo' + endfunction + + function foo = minmaxCEQ(z) + [tmp1,foo] = minmaxNonlinfun(z) + foo = foo' + endfunction + + dc = numderivative(minmaxC,z) + dceq = numderivative(minmaxCEQ,z) + endfunction + + // function handle defining new minmaxNonlinfun function + function [nc,nceq] = newNonlinfun(z) + [nc,nceq] = minmaxNonlinfun(z) + // add inequalities of the form Fi(x) - y <= 0 + tmp = [minmaxObjfun(z) - z(minmaxNumvar+1)]' + nc = [nc, tmp] + endfunction + + // function handle defining new gradient function for non-linear constraints + // this function passed when the gradient feature is on + function [dnc,dnceq] = newCGrad(z) + + // if constraint gradient present use it + if(flag2 == 1) then + [dnc, dnceq] = minmaxCGrad(z) + dnc = [dnc, zeros(size(dnc,'r'),1)] + dnceq = [dnceq, zeros(size(dnceq,'r'),1)] + else + // else use numderivative method to calculate gradient of constraints + [dnc, dnceq] = minmaxNonlinDer(z) + end + + // if objective gradient is present use it + if(flag1 == 1) then + derObjfun = minmaxFGrad(z) + mderObjfun = [derObjfun, -1*ones(size(derObjfun,'r'),1)] + dnc = [dnc; mderObjfun] + else + // else use numderivative to calculate gradient of set of obj functions + derObjfun = minmaxObjDer(z) + dnc = [dnc; derObjfun] + end + endfunction + + // to be passed as minimaxOptions to fmincon + if(flag1 == 1 | flag2 == 1) then + minmaxPassOptions = list("MaxIter", minmaxMaxIter, "CpuTime", minmaxCPU, "GradCon", newCGrad) + [x,fval,exitflag,output,lambda] = ... + fmincon(newObjfun,minmaxStartpoint,minmaxA,minmaxB,minmaxAeq,minmaxBeq,minmaxLb,minmaxUb,newNonlinfun,minmaxPassOptions) + + x = x(1:minmaxNumvar) + fval = minmaxObjfun(x) + maxfval = max(fval) + else + minmaxPassOptions = list("MaxIter", minmaxMaxIter, "CpuTime", minmaxCPU) + [x,fval,exitflag,output,lambda] = ... + fmincon(newObjfun,minmaxStartpoint,minmaxA,minmaxB,minmaxAeq,minmaxBeq,minmaxLb,minmaxUb,newNonlinfun,minmaxPassOptions) + + x = x(1:minmaxNumvar) + fval = minmaxObjfun(x) + maxfval = max(fval) + end + + +endfunction diff --git a/macros/fminimaxCheckdims.bin b/macros/fminimaxCheckdims.bin Binary files differnew file mode 100644 index 0000000..be1ed46 --- /dev/null +++ b/macros/fminimaxCheckdims.bin diff --git a/macros/fminimaxCheckdims.sci b/macros/fminimaxCheckdims.sci new file mode 100644 index 0000000..83cb5b5 --- /dev/null +++ b/macros/fminimaxCheckdims.sci @@ -0,0 +1,56 @@ + +// Copyright (C) 2010 - DIGITEO - Michael Baudin +// +// This file must be used under the terms of the GNU LGPL license. + +function errmsg = fminimaxCheckdims ( funname , var , varname , ivar , matdims ) + // Generates an error if the variable has not the required size. + // + // Calling Sequence + // errmsg = fminimaxCheckdims ( funname , var , varname , ivar , matdims ) + // + // Parameters + // funname : a 1 x 1 matrix of strings, the name of the calling function. + // var : a 1 x 1 matrix of valid Scilab data type, the variable + // varname : a 1 x 1 matrix of string, the name of the variable + // ivar : a 1 x 1 matrix of floating point integers, the index of the input argument in the calling sequence + // matdims : 1 x 2 matrix of floating point integers, the number of rows, columns for the variable #ivar + // errmsg : a 1 x 1 matrix of strings, the error message. If there was no error, the error message is the empty matrix. + // + // Description + // This function is designed to be used to design functions where + // the input argument has a known shape. + // This function cannot be use when var is a function, or more + // generally, for any input argument for which the size function + // does not work. + // Last update : 05/08/2010. + // + // Examples + // // The function takes a 2 x 3 matrix of doubles. + // function y = myfunction ( x ) + // fminimaxCheckdims ( "myfunction" , x , "x" , 1 , [2 3] ) + // y = x + // endfunction + // // Calling sequences which work + // y = myfunction ( ones(2,3) ) + // y = myfunction ( zeros(2,3) ) + // // Calling sequences which generate an error + // y = myfunction ( ones(1,3) ) + // y = myfunction ( zeros(2,4) ) + // + // Authors + // Michael Baudin - 2010 - DIGITEO + // + + [lhs,rhs]=argn() + fminimaxCheckrhs ( "fminimaxCheckdims" , rhs , 5 ) + fminimaxChecklhs ( "fminimaxCheckdims" , lhs , [0 1] ) + + errmsg = [] + if ( or ( size(var) <> matdims ) ) then + strexp = strcat(string(matdims)," ") + strcomp = strcat(string(size(var))," ") + errmsg = msprintf(gettext("%s: Expected size [%s] for input argument %s at input #%d, but got [%s] instead."), funname, strexp, varname , ivar , strcomp ); + error(errmsg) + end +endfunction diff --git a/macros/fminimaxChecklhs.bin b/macros/fminimaxChecklhs.bin Binary files differnew file mode 100644 index 0000000..7889084 --- /dev/null +++ b/macros/fminimaxChecklhs.bin diff --git a/macros/fminimaxChecklhs.sci b/macros/fminimaxChecklhs.sci new file mode 100644 index 0000000..78e4013 --- /dev/null +++ b/macros/fminimaxChecklhs.sci @@ -0,0 +1,79 @@ +// Copyright (C) 2010 - DIGITEO - Michael Baudin +// +// This file must be used under the terms of the GNU LGPL license. + +function errmsg = fminimaxChecklhs ( funname , lhs , lhsset ) + // Generates an error if the number of LHS is not in given set. + // + // Calling Sequence + // errmsg = fminimaxChecklhs ( funname , lhs , lhsset ) + // + // Parameters + // funname : a 1 x 1 matrix of strings, the name of the calling function. + // lhs : a 1 x 1 matrix of floating point integers, the actual number of output arguments + // lhsset : a 1 x n or n x 1 matrix of floating point integers, the authorized number of output arguments + // errmsg : a 1 x 1 matrix of strings, the error message. If there was no error, the error message is the empty matrix. + // + // Description + // This function is designed to be used to design functions with + // variable number of output arguments. + // Notice that it is useless to call this function if the + // function definition does not use the varargout statement. + // Notice that a function as a minimum of 1 output argument. + // Last update : 29/07/2010. + // + // Examples + // // The function takes 3 input arguments and 1/2 output arguments + // function varargout = myfunction ( x1 , x2 , x3 ) + // [lhs, rhs] = argn() + // fminimaxCheckrhs ( "myfunction" , rhs , 3 : 3 ) + // fminimaxChecklhs ( "myfunction" , lhs , 1 : 2 ) + // y1 = x1 + x2 + // y2 = x2 + x3 + // varargout(1) = y1 + // if ( lhs == 2 ) then + // varargout(2) = y2 + // end + // endfunction + // // Calling sequences which work + // myfunction ( 1 , 2 , 3 ) + // y1 = myfunction ( 1 , 2 , 3 ) + // [ y1 , y2 ] = myfunction ( 1 , 2 , 3 ) + // // Calling sequences which generate an error + // [ y1 , y2 , y3 ] = myfunction ( 1 , 2 , 3 ) + // + // // The function takes 1 or 3 output arguments, but not 2 + // function varargout = myfunction ( x1 , x2 , x3 ) + // [lhs, rhs] = argn() + // fminimaxCheckrhs ( "myfunction" , rhs , 3 : 3 ) + // fminimaxChecklhs ( "myfunction" , lhs , [1 3] ) + // y1 = x1 + x2 + // y2 = x2 + x3 + // y3 = x1 + x3 + // varargout(1) = y1 + // if ( lhs == 3 ) then + // varargout(2) = y2 + // varargout(3) = y3 + // end + // endfunction + // // Calling sequences which work + // myfunction ( 1 , 2 , 3 ) + // y1 = myfunction ( 1 , 2 , 3 ) + // [ y1 , y2 , y3 ] = myfunction ( 1 , 2 , 3 ) + // // Calling sequences which generate an error + // [y1 , y2] = myfunction ( 1 , 2 , 3 ) + // + // Authors + // Michael Baudin - 2010 - DIGITEO + // + + errmsg = [] + if ( and ( lhs <> lhsset ) ) then + lhsstr = strcat(string(lhsset)," ") + errmsg = msprintf(gettext("%s: Unexpected number of output arguments : %d provided while the expected number of output arguments should be in the set [%s]."), funname , lhs , lhsstr ); + error(errmsg) + end +endfunction + + + diff --git a/macros/fminimaxCheckrhs.bin b/macros/fminimaxCheckrhs.bin Binary files differnew file mode 100644 index 0000000..6b190f7 --- /dev/null +++ b/macros/fminimaxCheckrhs.bin diff --git a/macros/fminimaxCheckrhs.sci b/macros/fminimaxCheckrhs.sci new file mode 100644 index 0000000..7877153 --- /dev/null +++ b/macros/fminimaxCheckrhs.sci @@ -0,0 +1,102 @@ +// Copyright (C) 2010 - DIGITEO - Michael Baudin +// +// This file must be used under the terms of the GNU LGPL license. + +function errmsg = fminimaxCheckrhs ( funname , rhs , rhsset ) + // Generates an error if the number of RHS is not in given set. + // + // Calling Sequence + // errmsg = fminimaxCheckrhs ( funname , rhs , rhsset ) + // + // Parameters + // funname : a 1 x 1 matrix of strings, the name of the calling function. + // rhs : a 1 x 1 matrix of floating point integers, the actual number of input arguments + // rhsset : a 1 x n or n x 1 matrix of floating point integers, the authorized number of input arguments + // errmsg : a 1 x 1 matrix of strings, the error message. If there was no error, the error message is the empty matrix. + // + // Description + // This function is designed to be used to design functions with + // variable number of input arguments. + // Notice that it is useless to call this function if the + // function definition does not use the varargin statement. + // Last update : 05/08/2010. + // Last update : 29/07/2010. + // + // Examples + // // The function takes 2/3 input arguments and 1 output arguments + // function y = myfunction ( varargin ) + // [lhs, rhs] = argn() + // fminimaxCheckrhs ( "myfunction" , rhs , 2:3 ) + // fminimaxChecklhs ( "myfunction" , lhs , 1 ) + // x1 = varargin(1) + // x2 = varargin(2) + // if ( rhs >= 3 ) then + // x3 = varargin(3) + // else + // x3 = 2 + // end + // y = x1 + x2 + x3 + // endfunction + // // Calling sequences which work + // y = myfunction ( 1 , 2 ) + // y = myfunction ( 1 , 2 , 3 ) + // // Calling sequences which generate an error + // y = myfunction ( 1 ) + // y = myfunction ( 1 , 2 , 3 , 4 ) + // + // // The function takes 2 or 4 input arguments, but not 3 + // function y = myfunction ( varargin ) + // [lhs, rhs] = argn() + // fminimaxCheckrhs ( "myfunction" , rhs , [2 4] ) + // fminimaxChecklhs ( "myfunction" , lhs , 1 ) + // x1 = varargin(1) + // x2 = varargin(2) + // if ( rhs >= 3 ) then + // x3 = varargin(3) + // x4 = varargin(4) + // else + // x3 = 2 + // x4 = 3 + // end + // y = x1 + x2 + x3 + x4 + // endfunction + // // Calling sequences which work + // y = myfunction ( 1 , 2 ) + // y = myfunction ( 1 , 2 , 3 , 4 ) + // // Calling sequences which generate an error + // y = myfunction ( 1 ) + // y = myfunction ( 1 , 2 , 3 ) + // y = myfunction ( 1 , 2 , 3 , 4, 5 ) + // + // // The function takes 2 input arguments and 0/1 output arguments. + // // Notice that if the checkrhs function is not called, + // // the variable x2 might be used from the user's context, + // // that is, if the caller has defined the variable x2, it + // // is used in the callee. + // // Here, we want to avoid this. + // function y = myfunction ( x1 , x2 ) + // [lhs, rhs] = argn() + // fminimaxCheckrhs ( "myfunction" , rhs , 2 ) + // fminimaxChecklhs ( "myfunction" , lhs , [0 1] ) + // y = x1 + x2 + // endfunction + // // Calling sequences which work + // y = myfunction ( 1 , 2 ) + // // Calling sequences which generate an error + // y = myfunction ( 1 ) + // y = myfunction ( 1 , 2 , 3 ) + // + // Authors + // Michael Baudin - 2010 - DIGITEO + // + + errmsg = [] + if ( and(rhs <> rhsset) ) then + rhsstr = strcat(string(rhsset)," ") + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while the number of expected input arguments should be in the set [%s]."), funname , rhs , rhsstr ); + error(errmsg) + end +endfunction + + + diff --git a/macros/fminimaxChecktype.bin b/macros/fminimaxChecktype.bin Binary files differnew file mode 100644 index 0000000..4f42fe0 --- /dev/null +++ b/macros/fminimaxChecktype.bin diff --git a/macros/fminimaxChecktype.sci b/macros/fminimaxChecktype.sci new file mode 100644 index 0000000..48514b6 --- /dev/null +++ b/macros/fminimaxChecktype.sci @@ -0,0 +1,65 @@ +// Copyright (C) 2010 - DIGITEO - Michael Baudin +// +// This file must be used under the terms of the GNU LGPL license. + +function errmsg = fminimaxChecktype ( funname , var , varname , ivar , expectedtype ) + // Generates an error if the given variable is not of expected type. + // + // Calling Sequence + // errmsg = fminimaxChecktype ( funname , var , varname , ivar , expectedtype ) + // + // Parameters + // funname : a 1 x 1 matrix of strings, the name of the calling function. + // var : a 1 x 1 matrix of valid Scilab data type, the variable + // varname : a 1 x 1 matrix of string, the name of the variable + // ivar : a 1 x 1 matrix of floating point integers, the index of the input argument in the calling sequence + // expectedtype : a n x 1 or 1 x n matrix of strings, the available types for the variable #ivar + // errmsg : a 1 x 1 matrix of strings, the error message. If there was no error, the error message is the empty matrix. + // + // Description + // This function is designed to be used to design functions with + // input arguments with variable type. + // We use the typeof function to compute the type of the variable: + // see help typeof to get the list of all available values for expectedtype. + // Last update : 29/07/2010. + // + // Examples + // // The function takes a string argument. + // function myfunction ( x ) + // fminimaxChecktype ( "myfunction" , x , "x" , 1 , "string" ) + // disp("This is a string") + // endfunction + // // Calling sequences which work + // myfunction ( "Scilab" ) + // // Calling sequences which generate an error + // myfunction ( 123456 ) + // + // // The function takes a string or a matrix of doubles argument. + // function myfunction ( x ) + // fminimaxChecktype ( "myfunction" , x , "x" , 1 , [ "string" "constant" ] ) + // if ( typeof(x) == "string" ) then + // disp("This is a matrix of strings") + // else + // disp("This is a matrix of doubles") + // end + // endfunction + // // Calling sequences which work + // myfunction ( "Scilab" ) + // myfunction ( 123456 ) + // // Calling sequences which generate an error + // myfunction ( uint8(2) ) + // + // Authors + // Michael Baudin - 2010 - DIGITEO + // + + errmsg = [] + if ( and ( typeof ( var ) <> expectedtype ) ) then + strexp = """" + strcat(expectedtype,""" or """) + """" + errmsg = msprintf(gettext("%s: Expected type [%s] for input argument %s at input #%d, but got ""%s"" instead."),funname, strexp, varname , ivar , typeof(var) ); + error(errmsg); + end +endfunction + + + diff --git a/macros/fminimaxCheckvector.bin b/macros/fminimaxCheckvector.bin Binary files differnew file mode 100644 index 0000000..28af2c5 --- /dev/null +++ b/macros/fminimaxCheckvector.bin diff --git a/macros/fminimaxCheckvector.sci b/macros/fminimaxCheckvector.sci new file mode 100644 index 0000000..87c8553 --- /dev/null +++ b/macros/fminimaxCheckvector.sci @@ -0,0 +1,63 @@ +// Copyright (C) 2010 - DIGITEO - Michael Baudin +// +// This file must be used under the terms of the GNU LGPL license. + +function errmsg = fminimaxCheckvector ( funname , var , varname , ivar , nbval ) + // Generates an error if the variable is not a vector. + // + // Calling Sequence + // errmsg = fminimaxCheckvector ( funname , var , varname , ivar ) + // + // Parameters + // funname : a 1 x 1 matrix of strings, the name of the calling function. + // var : a 1 x 1 matrix of valid Scilab data type, the variable + // varname : a 1 x 1 matrix of string, the name of the variable + // ivar : a 1 x 1 matrix of floating point integers, the index of the input argument in the calling sequence + // nbval : a 1 x 1 matrix of floating point integers, the number of entries in the vector. + // errmsg : a 1 x 1 matrix of strings, the error message. If there was no error, the error message is the empty matrix. + // + // Description + // This function is designed to be used to design functions where + // the input argument is a vector, that is, a matrix for which + // nrows == 1 or ncols == 1. + // This function cannot be use when var is a function, or more + // generally, for any input argument for which the size function + // does not work. + // + // Examples + // // The function takes a vector of 3 doubles. + // function y = myfunction ( x ) + // fminimaxCheckvector ( "myfunction" , x , "x" , 1 , 3 ) + // y = x + // endfunction + // // Calling sequences which work + // y = myfunction ( ones(1,3) ) + // y = myfunction ( zeros(3,1) ) + // // Calling sequences which generate an error + // // The following are not vectors + // y = myfunction ( ones(2,3) ) + // y = myfunction ( zeros(3,2) ) + // // The following have the wrong number of entries + // y = myfunction ( ones(1,3) ) + // + // Authors + // Michael Baudin - 2010 - DIGITEO + // + + errmsg = [] + nrows = size(var,"r") + ncols = size(var,"c") + if ( nrows <> 1 & ncols <> 1 ) then + strcomp = strcat(string(size(var))," ") + errmsg = msprintf(gettext("%s: Expected a vector matrix for input argument %s at input #%d, but got [%s] instead."), funname, varname , ivar , strcomp ); + error(errmsg) + end + if ( ( nrows == 1 & ncols <> nbval ) | ( ncols == 1 & nrows <> nbval ) ) then + strcomp = strcat(string(size(var))," ") + errmsg = msprintf(gettext("%s: Expected %d entries for input argument %s at input #%d, but current dimensions are [%s] instead."), funname, nbval , varname , ivar , strcomp ); + error(errmsg) + end +endfunction + + + diff --git a/macros/fminunc.bin b/macros/fminunc.bin Binary files differnew file mode 100644 index 0000000..b65c851 --- /dev/null +++ b/macros/fminunc.bin 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 Binary files differdiff --git a/macros/linprog.bin b/macros/linprog.bin Binary files differnew file mode 100644 index 0000000..b9ba12a --- /dev/null +++ b/macros/linprog.bin diff --git a/macros/linprog.sci b/macros/linprog.sci new file mode 100644 index 0000000..4c949ba --- /dev/null +++ b/macros/linprog.sci @@ -0,0 +1,199 @@ +// 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: Guru Pradeep Reddy, Bhanu Priya Sayal +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in + + +function [xopt,fopt,exitflag,output,lambda] = linprog (varargin) + // Solves a linear programming problem. + // + // Calling Sequence + // xopt = linprog(c,A,b) + // xopt = linprog(c,A,b,Aeq,beq) + // xopt = linprog(c,A,b,Aeq,beq,lb,ub) + // xopt = linprog(c,A,b,Aeq,beq,lb,ub,param) + // [xopt, fopt, exitflag, output, lambda] = linprog(file) + // [xopt,fopt,exitflag,output,lambda] = linprog( ... ) + // + // Parameters + // c : a vector of double, contains coefficients of the variables in the objective + // A : a matrix of double, represents the linear coefficients in the inequality constraints A⋅x ≤ b. + // b : a vector of double, represents the linear coefficients in the inequality constraints A⋅x ≤ b. + // Aeq : a matrix of double, represents the linear coefficients in the equality constraints Aeq⋅x = beq. + // beq : a vector of double, represents the linear coefficients in the equality constraints Aeq⋅x = beq. + // lb : Lower bounds, specified as a vector or array of double. lb represents the lower bounds elementwise in lb ≤ x ≤ ub. + // ub : Upper bounds, specified as a vector or array of double. ub represents the upper bounds elementwise in lb ≤ x ≤ ub. + // options : a list containing the parameters to be set. + // file : a string describing the path to the mps file. + // xopt : a vector of double, the computed solution of the optimization problem. + // fopt : a double, the value of the function at x. + // status : status flag returned from symphony. See below for details. + // output : The output data structure contains detailed information about the optimization process. See below for details. + // lambda : The structure consist of the Lagrange multipliers at the solution of problem. See below for details. + // + // Description + // OSI-CLP is used for solving the linear programming problems, OSI-CLP is a library written in C++. + // Search the minimum of a constrained linear programming problem specified by : + // + // <latex> + // \begin{eqnarray} + // &\mbox{min}_{x} + // & c^T⋅x \\ + // & \text{subject to} & A⋅x \leq b \\ + // & & Aeq⋅x = beq \\ + // & & lb \leq x \leq ub \\ + // \end{eqnarray} + // </latex> + // The routine calls Clp for solving the linear programming problem, Clp is a library written in C++. + // + // 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 : Primal Infeasible </listitem> + // <listitem>exitflag=2 : Dual Infeasible</listitem> + // <listitem>exitflag=3 : Maximum Number of Iterations Exceeded. Output may not be optimal.</listitem> + // <listitem>exitflag=4 : Solution Abandoned</listitem> + // <listitem>exitflag=5 : Primal objective limit reached.</listitem> + // <listitem>exitflag=6 : Dual objective limit reached.</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.constrviolation: The max-norm of the constraint violation.</listitem> + // </itemizedlist> + // + // The lambda data structure contains the Lagrange multipliers at the end + // of optimization. In the current version the values are returned only when the the solution is optimal. + // It has type "struct" and contains the following fields. + // <itemizedlist> + // <listitem>lambda.lower: The Lagrange multipliers for the lower bound constraints.</listitem> + // <listitem>lambda.upper: The Lagrange multipliers for the upper bound constraints.</listitem> + // <listitem>lambda.eqlin: The Lagrange multipliers for the linear equality constraints.</listitem> + // <listitem>lambda.ineqlin: The Lagrange multipliers for the linear inequality constraints.</listitem> + // </itemizedlist> + // + // Examples + // //Optimal problems + // //Linear program, linear inequality constraints + // c=[-1,-1/3]' + // A=[1,1;1,1/4;1,-1;-1/4,-1;-1,-1;-1,1] + // b=[2,1,2,1,-1,2] + // [xopt,fopt,exitflag,output,lambda]=linprog(c, A, b) + // // Press ENTER to continue + // + // Examples + // //Linear program with Linear Inequalities and Equalities` + // c=[-1,-1/3]' + // A=[1,1;1,1/4;1,-1;-1/4,-1;-1,-1;-1,1] + // b=[2,1,2,1,-1,2] + // Aeq=[1,1/4] + // beq=[1/2] + // [xopt,fopt,exitflag,output,lambda]=linprog(c, A, b, Aeq, beq) + // // Press ENTER to continue + // + // Examples + // //Linear program with all constraint types + // c=[-1,-1/3]' + // A=[1,1;1,1/4;1,-1;-1/4,-1;-1,-1;-1,1] + // b=[2,1,2,1,-1,2] + // Aeq=[1,1/4] + // beq=[1/2] + // lb=[-1,-0.5] + // ub=[1.5,1.25] + // [xopt,fopt,exitflag,output,lambda]=linprog(c, A, b, Aeq, beq, lb, ub) + // // Press ENTER to continue + // + // Examples + // //Primal Infeasible Problem + // c=[-1,-1,-1]' + // A=[1,2,-1] + // b=[-4] + // Aeq=[1,5,3;1,1,0] + // beq=[10,100] + // lb=[0,0,0] + // ub=[%inf,%inf,%inf] + // [xopt,fopt,exitflag,output,lambda]= linprog(c,A,b,Aeq,beq,lb,ub) + // // Press ENTER to continue + // + // Examples + // //Dual Infeasible Problem + // c=[3,5,-7]' + // A=[-1,-1,4;1,1,4] + // b=[-8,5] + // Aeq=[] + // beq=[] + // lb=[-%inf,-%inf,-%inf] + // ub=[%inf,%inf,%inf] + // [xopt,fopt,exitflag,output,lambda]= linprog(c,A,b,Aeq,beq,lb,ub) + // // Press ENTER to continue + // + // Examples + // filepath = get_absolute_file_path('linprog.dem.sce'); + // filepath = filepath + "exmip1.mps" + // [xopt,fopt,exitflag,output,lambda] =linprog(filepath) + // Authors + // Bhanu Priya Sayal, Guru Pradeep Reddy + + if(type(varargin(1))==1) then + [lhs , rhs] = argn(); + //To check the number of argument given by user + if ( rhs < 3 | rhs == 4 | rhs == 6 | rhs >8 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be in the set of [3 5 7 8]"), "linprog", rhs); + error(errmsg) + end + + c = varargin(1); + A = varargin(2); + b = varargin(3); + + if ( rhs<4 ) then + Aeq = [] + beq = [] + else + Aeq = varargin(4); + beq = varargin(5); + end + + if ( rhs<6 ) then + lb = []; + ub = []; + else + lb = varargin(6); + ub = varargin(7); + end + + if ( rhs<8 | size(varargin(8)) ==0 ) then + param = list(); + else + param =varargin(8); + end + [xopt,fopt,exitflag,output,lambda]=matrix_linprog(c,A,b,Aeq,beq,lb,ub,param); + elseif(type(varargin(1))==10) then + + [lhs , rhs] = argn(); + + //To check the number of argument given by user + if ( rhs < 1 | rhs > 2) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be in the set of [1 2]"),"linprog",rhs); + error(errmsg) + end + mpsFile = varargin(1); + if ( rhs<2 | size(varargin(2)) ==0 ) then + param = list(); + else + param =varargin(2); + end + [xopt,fopt,exitflag,output,lambda]=mps_linprog(mpsFile,param); + end + +endfunction diff --git a/macros/lsqlin.bin b/macros/lsqlin.bin Binary files differindex eb02534..7018f4a 100644 --- a/macros/lsqlin.bin +++ b/macros/lsqlin.bin diff --git a/macros/lsqlin.sci b/macros/lsqlin.sci index 960b4db..97aaae6 100644 --- a/macros/lsqlin.sci +++ b/macros/lsqlin.sci @@ -151,6 +151,11 @@ function [xopt,resnorm,residual,exitflag,output,lambda] = lsqlin (varargin) b = varargin(4); nbVar = size(C,2); + if(nbVar == 0) then + errmsg = msprintf(gettext("%s: Cannot determine the number of variables because input objective coefficients is empty"), "lsqlin"); + error(errmsg); + end + if ( rhs<5 ) then Aeq = [] beq = [] diff --git a/macros/lsqnonneg.bin b/macros/lsqnonneg.bin Binary files differindex 7ea7f4a..39c553b 100644 --- a/macros/lsqnonneg.bin +++ b/macros/lsqnonneg.bin diff --git a/macros/lsqnonneg.sci b/macros/lsqnonneg.sci index 3b152bf..3e5ea81 100644 --- a/macros/lsqnonneg.sci +++ b/macros/lsqnonneg.sci @@ -97,6 +97,12 @@ function [xopt,resnorm,residual,exitflag,output,lambda] = lsqnonneg (varargin) C = varargin(1); d = varargin(2); nbVar = size(C,2); + + if(nbVar == 0) then + errmsg = msprintf(gettext("%s: Cannot determine the number of variables because input objective coefficients is empty"), "lsqnonneg"); + error(errmsg); + end + if ( rhs<3 | size(varargin(3)) ==0 ) then param = list(); else diff --git a/macros/matrix_linprog.bin b/macros/matrix_linprog.bin Binary files differnew file mode 100644 index 0000000..21e525c --- /dev/null +++ b/macros/matrix_linprog.bin diff --git a/macros/matrix_linprog.sci b/macros/matrix_linprog.sci new file mode 100755 index 0000000..1f6a6bf --- /dev/null +++ b/macros/matrix_linprog.sci @@ -0,0 +1,245 @@ +// Copyright (C) 2015 - IIT Bombay - FOSSEE +// +// Author: Guru Pradeep Reddy, Bhanu Priya Sayal +// Organization: FOSSEE, IIT Bombay +// Email: gurupradeept@gmail.com, bhanupriyasayal@gmail.com +// 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 + + +function [xopt,fopt,exitflag,output,lambda] = matrix_linprog (varargin) +//To check the number of input and output argument + [lhs , rhs] = argn(); + + c = []; + A = []; + b = []; + Aeq = []; + beq = []; + lb = []; + ub = []; + options = list(); + +//To check the number of argument given by user + if ( rhs < 3 | rhs == 4 | rhs == 6 | rhs >8 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be in the set of [3 5 7 8]"), "linprog", rhs); + error(errmsg) + end + + c = varargin(1); + + if (size(c,2)~=1) then + errmsg = msprintf(gettext("%s: Objective Coefficients should be a column matrix"), "linprog"); + error(errmsg); + end + + if(size(c,2) == 0) then + errmsg = msprintf(gettext("%s: Cannot determine the number of variables because input objective coefficients is empty"), "linprog"); + error(errmsg); + end + + nbVar = size(c,1); + A = varargin(2); + b = varargin(3); + + if ( rhs<4 ) then + Aeq = [] + beq = [] + else + Aeq = varargin(4); + beq = varargin(5); + end + + if ( rhs<6 ) then + lb = []; + ub = []; + else + lb = varargin(6); + ub = varargin(7); + end + + if ( rhs<8 | size(varargin(8)) ==0 ) then + param = list(); + else + param =varargin(8); + end + + nbConInEq = size(A,1); + nbConEq = size(Aeq,1); + + if (size(lb,2)== [nbVar]) then + lb = lb'; + end + + if (size(ub,2)== [nbVar]) then + ub = ub'; + end + + if (size(b,2)== [nbConInEq]) then + b = b'; + end + + if (size(beq,2)== [nbConEq]) then + beq = beq'; + end + + + if (size(lb,2)==0) then + lb = repmat(-%inf,nbVar,1); + end + + if (size(ub,2)==0) then + ub = repmat(%inf,nbVar,1); + end + + if (type(param) ~= 15) then + errmsg = msprintf(gettext("%s: options should be a list "), "linprog"); + error(errmsg); + end + + + if (modulo(size(param),2)) then + errmsg = msprintf(gettext("%s: Size of parameters should be even"), "linprog"); + error(errmsg); + end + + options = list("MaxIter", [3000]); + + for i = 1:(size(param))/2 + + select param(2*i-1) + case "MaxIter" then + options(2*i) = param(2*i); + else + errmsg = msprintf(gettext("%s: Unrecognized parameter name ''%s''."), "qpipoptmat", param(2*i-1)); + error(errmsg) + end + end + + //Check the size of inequality constraint which should be equal to the number of variables + if ( size(A,2) ~= nbVar & size(A,2) ~= 0) then + errmsg = msprintf(gettext("%s: The number of columns in A must be the same as the number of elements of c"), "linprog"); + error(errmsg); + end + + //Check the size of equality constraint which should be equal to the number of variables + if ( size(Aeq,2) ~= nbVar & size(Aeq,2) ~= 0 ) then + errmsg = msprintf(gettext("%s: The number of columns in Aeq must be the same as the number of elements of c"), "linprog"); + error(errmsg); + end + + //Check the size of Lower Bound which should be equal to the number of variables + if ( size(lb,1) ~= nbVar) then + errmsg = msprintf(gettext("%s: The Lower Bound is not equal to the number of variables"), "linprog"); + error(errmsg); + end + + //Check the size of Upper Bound which should equal to the number of variables + if ( size(ub,1) ~= nbVar) then + errmsg = msprintf(gettext("%s: The Upper Bound is not equal to the number of variables"), "linprog"); + error(errmsg); + end + //Check the size of constraints of Lower Bound which should equal to the number of constraints + if ( size(b,1) ~= nbConInEq & size(b,2) ~= 0) then + errmsg = msprintf(gettext("%s: The number of rows in A must be the same as the number of elements of b"), "linprog"); + error(errmsg); + end + + //Check the size of constraints of Upper Bound which should equal to the number of constraints + if ( size(beq,1) ~= nbConEq & size(beq,2) ~= 0) then + errmsg = msprintf(gettext("%s: The number of rows in Aeq must be the same as the number of elements of beq"), "linprog"); + error(errmsg); + end + + //Check if the user gives a matrix instead of a vector + + if (size(lb,1)~=1)& (size(lb,2)~=1) then + errmsg = msprintf(gettext("%s: Lower Bound should be a vector"), "linprog"); + error(errmsg); + end + + if (size(ub,1)~=1)& (size(ub,2)~=1) then + errmsg = msprintf(gettext("%s: Upper Bound should be a vector"), "linprog"); + error(errmsg); + end + + if (nbConInEq) then + if ((size(b,1)~=1)& (size(b,2)~=1)) then + errmsg = msprintf(gettext("%s: Constraint Lower Bound should be a vector"), "linprog"); + error(errmsg); + end + end + + if (nbConEq) then + if (size(beq,1)~=1)& (size(beq,2)~=1) then + errmsg = msprintf(gettext("%s: Constraint should be a vector"), "linprog"); + error(errmsg); + end + end + + for i = 1:nbConInEq + if (b(i) == -%inf) + errmsg = msprintf(gettext("%s: Value of b can not be negative infinity"), "linprog"); + error(errmsg); + end + end + + for i = 1:nbConEq + if (beq(i) == -%inf) + errmsg = msprintf(gettext("%s: Value of beq can not be negative infinity"), "linprog"); + error(errmsg); + end + end + + lb = lb(:); + ub = ub(:) + b = b(:); + beq = beq(:); + nbVar = size(c,1); + c = c'; + conMatrix = [Aeq;A]; + nbCon = size(conMatrix,1); + conlb = [beq; repmat(-%inf,nbConInEq,1)]; + conub = [beq;b]; + + [xopt,fopt,status,iter,Zl,dual] = linearprog(nbVar,nbCon,c,conMatrix,conlb,conub,lb',ub',options); + + xopt = xopt'; + exitflag = status; + output = struct("Iterations" , [],.. + "constrviolation" , []); + + output.Iterations = iter; + output.constrviolation = max([0;norm(Aeq*xopt-beq, 'inf');(lb-xopt);(xopt-ub);(A*xopt-b)]); + lambda = struct("reduced_cost" , [], .. + "ineqlin" , [], .. + "eqlin" , []); + + lambda.reduced_cost = Zl; + lambda.eqlin = dual(1:nbConEq); + lambda.ineqlin = dual(nbConEq+1:nbCon); + select status + + case 0 then + printf("\nOptimal Solution.\n"); + case 1 then + printf("\nPrimal Infeasible.\n"); + case 2 then + printf("\nDual Infeasible.\n"); + case 3 then + printf("\nIteration limit reached.\n"); + case 4 then + printf("\nNumerical Difficulties.\n"); + case 5 then + printf("\nPrimal Objective limit reached.\n"); + case 6 then + printf("\nDual Objective limit reached.\n"); + else + printf("\nInvalid status returned. Notify the Toolbox authors\n"); + break; + end + +endfunction diff --git a/macros/mps_linprog.bin b/macros/mps_linprog.bin Binary files differnew file mode 100644 index 0000000..ee7d822 --- /dev/null +++ b/macros/mps_linprog.bin diff --git a/macros/mps_linprog.sci b/macros/mps_linprog.sci new file mode 100644 index 0000000..55f3edf --- /dev/null +++ b/macros/mps_linprog.sci @@ -0,0 +1,90 @@ +// Copyright (C) 2015 - IIT Bombay - FOSSEE +// +// Author: Bhanu Priya Sayal, Guru Pradeep Reddy +// Organization: FOSSEE, IIT Bombay +// Email:bhanupriyasayal@gmail.com,gurupradeept@gmail.com +// 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 + + +function [xopt,fopt,exitflag,output,lambda] =mps_linprog(varargin) + + //To check the number of input and output argument + [lhs , rhs] = argn(); + + //To check the number of argument given by user + if ( rhs < 1 | rhs > 2) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be in the set of [1 2]"),"linprog",rhs); + error(errmsg) + end + mpsFile = varargin(1); + + [sizeFile isFile] = fileinfo(mpsFile); + if(isFile ~= 0) then + errmsg = msprintf(gettext("%s: File is not present at the given location"),"linprog",rhs); + error(errmsg) + end + + if ( rhs<2 | size(varargin(2)) ==0 ) then + param = list(); + else + param =varargin(2); + end + + if (type(param) ~= 15) then + errmsg = msprintf(gettext("%s: options should be a list "), "mps_linprog"); + error(errmsg); + end + + if (modulo(size(param),2)) then + errmsg = msprintf(gettext("%s: Size of parameters should be even"), "mps_linprog"); + error(errmsg); + end + options = list("MaxIter" , [3000],); + + for i = 1:(size(param))/2 + select param(2*i-1) + case "MaxIter" then + options(2*i) = param(2*i); + end + end + //Calling the function by passing the required parameters + + [xopt,fopt,status,iter,Zl,dual] = rmps(mpsFile,options); + + xopt = xopt'; + fopt=fopt; + exitflag = status; + output = struct("Iterations" , []); + output.Iterations = iter; + lambda = struct("reduced_cost" , [], .. + "dual" ,[]); + + lambda.reduced_cost = Zl; + lambda.dual =dual; + + select status + + case 0 then + printf("\nOptimal Solution.\n"); + case 1 then + printf("\nPrimal Infeasible.\n"); + case 2 then + printf("\nDual Infeasible.\n"); + case 3 then + printf("\nIteration limit reached.\n"); + case 4 then + printf("\nNumerical Difficulties.\n"); + case 5 then + printf("\nPrimal Objective limit reached.\n"); + case 6 then + printf("\nDual Objective limit reached.\n"); + else + printf("\nInvalid status returned. Notify the Toolbox authors\n"); + break; + end + +endfunction diff --git a/macros/names b/macros/names index e04f783..1055c2a 100644 --- a/macros/names +++ b/macros/names @@ -1,5 +1,19 @@ +fgoalattain +fgoalattainFunctions +fminbnd +fmincon +fminimax +fminimaxCheckdims +fminimaxChecklhs +fminimaxCheckrhs +fminimaxChecktype +fminimaxCheckvector +fminunc +linprog lsqlin lsqnonneg +matrix_linprog +mps_linprog qpipopt qpipoptmat setOptions diff --git a/macros/qpipoptmat.bin b/macros/qpipoptmat.bin Binary files differindex 64e0b44..58424ec 100644 --- a/macros/qpipoptmat.bin +++ b/macros/qpipoptmat.bin diff --git a/macros/qpipoptmat.sci b/macros/qpipoptmat.sci index b74c718..ca32f8e 100644 --- a/macros/qpipoptmat.sci +++ b/macros/qpipoptmat.sci @@ -145,6 +145,10 @@ function [xopt,fopt,exitflag,output,lambda] = qpipoptmat (varargin) f = varargin(2); nbVar = size(H,1); + if(nbVar == 0) then + errmsg = msprintf(gettext("%s: Cannot determine the number of variables because input objective coefficients is empty"), "qpipoptmat"); + error(errmsg); + end if ( rhs<3 ) then A = [] @@ -424,5 +428,4 @@ function [xopt,fopt,exitflag,output,lambda] = qpipoptmat (varargin) break; end - endfunction diff --git a/macros/symphonymat.bin b/macros/symphonymat.bin Binary files differindex ec35ee2..a5e68ab 100644 --- a/macros/symphonymat.bin +++ b/macros/symphonymat.bin diff --git a/macros/symphonymat.sci b/macros/symphonymat.sci index 2d51b84..3848850 100644 --- a/macros/symphonymat.sci +++ b/macros/symphonymat.sci @@ -193,12 +193,16 @@ function [xopt,fopt,status,iter] = symphonymat (varargin) A = varargin(3) b = varargin(4) + if(size(c,2) == 0) then + errmsg = msprintf(gettext("%s: Cannot determine the number of variables because input objective coefficients is empty"),"Symphonymat"); + error(errmsg); + end + if (size(c,2)~=1) then errmsg = msprintf(gettext("%s: Objective Coefficients should be a column matrix"), "Symphonymat"); error(errmsg); end - nbVar = size(c,1); if ( rhs<5 ) then |