// 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
//
//
// \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}
//
//
// 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.
//
// Syntax : options= list("MaxIter", [---], "CpuTime", [---], "GradObj", ---, "Hessian", ---, "GradCon", ---);
// MaxIter : a Scalar, containing the Maximum Number of Iteration that the solver should take.
// CpuTime : a Scalar, containing the Maximum amount of CPU Time that the solver should take.
// GradObj : a function, representing the gradient function of the Objective in Vector Form.
// 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.
// 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.
// Default Values : options = list("MaxIter", [3000], "CpuTime", [600]);
//
//
// The exitflag allows to know the status of the optimization which is given back by Ipopt.
//
// exitflag=0 : Optimal Solution Found
// exitflag=1 : Maximum Number of Iterations Exceeded. Output may not be optimal.
// exitflag=2 : Maximum amount of CPU Time exceeded. Output may not be optimal.
// exitflag=3 : Stop at Tiny Step.
// exitflag=4 : Solved To Acceptable Level.
// exitflag=5 : Converged to a point of local infeasibility.
//
//
// 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.
//
// output.Iterations: The number of iterations performed during the search
// output.Cpu_Time: The total cpu-time spend during the search
// output.Objective_Evaluation: The number of Objective Evaluations performed during the search
// output.Dual_Infeasibility: The Dual Infeasiblity of the final soution
// output.Message: The output message for the problem
//
//
// 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.
//
// lambda.lower: The Lagrange multipliers for the lower bound constraints.
// lambda.upper: The Lagrange multipliers for the upper bound constraints.
// lambda.eqlin: The Lagrange multipliers for the linear equality constraints.
// lambda.ineqlin: The Lagrange multipliers for the linear inequality constraints.
// lambda.eqnonlin: The Lagrange multipliers for the non-linear equality constraints.
// lambda.ineqnonlin: The Lagrange multipliers for the non-linear inequality constraints.
//
//
// 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)
// // Press ENTER to continue
// 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>10 ) then
errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while it should be 4,6,8,9,10"), "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,10s"), "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 (fun) is a function or not
Checktype("fmincon", fun, "f", 1, "function");
//To check whether the 2nd Input argument (x0) is a vector/scalar
Checktype("fmincon", x0, "x0", 2, "constant");
//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
Checktype("fmincon", A, "A", 3, "constant");
//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
Checktype("fmincon", b, "b", 4, "constant");
//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
Checktype("fmincon", Aeq, "Aeq", 5, "constant");
//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
Checktype("fmincon", beq, "beq", 6, "constant");
//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
Checktype("fmincon", lb, "lb", 7, "constant");
//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
Checktype("fmincon", ub, "ub", 8, "constant");
//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'
Checktype("fmincon", param, "options", 10, "list");
//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
//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 convstr(param(2*i-1),'l')
case "maxiter" then
if (type(param(2*i))~=1) then
errmsg = msprintf(gettext("%s: Value for Maximum Iteration should be a Constant"), "fmincon");
error(errmsg);
else
options(2) = param(2*i); //Setting the maximum number of iterations as per user entry
end
case "cputime" then
if (type(param(2*i))~=1) then
errmsg = msprintf(gettext("%s: Value for Maximum Cpu-time should be a Constant"), "fmincon");
error(errmsg);
else
options(4) = param(2*i); //Setting the maximum CPU time as per user entry
end
case "gradobj" then
if (type(param(2*i))==10) then
if (convstr(param(2*i))=="off") then
flag1 =0;
else
errmsg = msprintf(gettext("%s: Unrecognized String [%s] entered for the option- %s."), "fmincon",param(2*i), param(2*i-1));
error(errmsg);
end
else
flag1 = 1;
fGrad = param(2*i);
end
case "hessian" then
if (type(param(2*i))==10) then
if (convstr(param(2*i))=="off") then
flag2 =0;
else
errmsg = msprintf(gettext("%s: Unrecognized String [%s] entered for the option- %s."), "fmincon",param(2*i), param(2*i-1));
error(errmsg);
end
else
flag2 = 1;
lHess = param(2*i);
end
case "gradcon" then
if (type(param(2*i))==10) then
if (convstr(param(2*i))=="off") then
flag3 =0;
else
errmsg = msprintf(gettext("%s: Unrecognized String [%s] entered for the option- %s."), "fmincon",param(2*i), param(2*i-1));
error(errmsg);
end
else
flag3 = 1;
cGrad = param(2*i);
end
else
errmsg = msprintf(gettext("%s: Unrecognized parameter name %s."), "fmincon", param(2*i-1));
error(errmsg);
end
end
//To check for correct input of Objective Gradient function from the user
if (flag1==1) then
Checktype("fmincon", fGrad, "fGrad", 10, "function");
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
else
errmsg = msprintf(gettext("%s: Wrong Input for Objective Gradient function(3rd Parameter)---->Row Vector function of size [1 X %d] is Expected"), "fmincon",s(2));
error(errmsg);
end
end
//To check for correct input of Lagrangian Hessian function from the user
if (flag2==1) then
Checktype("fmincon", lHess, "lHess", 10, "function");
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
end
//To check for correct input of Constraint Gradient function from the user
if (flag3==1) then
Checktype("fmincon", cGrad, "cGrad", 10, "function");
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
end
//Defining an inbuilt Objective gradient function
function [y,check] = fGrad1(x)
if flag1==1 then
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
else
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
end
endfunction
//Defining an inbuilt Lagrangian Hessian function
function [y,check] = lHess1(x,obj,lambda)
if flag2==1 then
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
else
if(execstr('[grad,y1]=numderivative(fun,x)','errcatch')==10000)
check1=1;
else
[grad,y1]=numderivative(fun,x);
if (isreal(y1)==%F) then
check1=1;
else
check1=0;
end
end
if check1==0 then
if no_nlc~=0 then
if(execstr('[grad,y2]=numderivative(addnlc,x)','errcatch')==10000)
check2=1;
else
[grad,y2]=numderivative(addnlc,x);
if (isreal(y2)==%F) then
check2=1;
else
check2=0;
end
end
if check2==0 then
y2=matrix(y2, no_nlc*s(2)*s(2),1)
for i = 1:s(2)*s(2)
y(i)=0;
for j = 1:no_nlc
y(i)= y(i) + lambda(j)*y2((i-1)*no_nlc+j);
end
end
for i = 1:s(2)*s(2)
y(i) = y(i)+ obj*y1(i);
end
check=0;
else
check=1;
end
else
check=0;
for i = 1:s(2)*s(2)
y(i) = obj*y1(i);
end
end
else
check=1;
y=[0];
end
end
endfunction
//Defining an inbuilt Constraint gradient function
function [y,check] = addcGrad1(x)
if flag3==1 then
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
else
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
end
endfunction
//Creating a Dummy Variable for IPopt use
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,A,b,Aeq,beq,lb,ub,no_nlc,no_nlic,addnlc1,fGrad1,lHess1,addcGrad1,x0,options,empty)
//Calculating the values for the output
xopt = xopt';
exitflag = status;
output = struct("Iterations", [],"Cpu_Time",[],"Objective_Evaluation",[],"Dual_Infeasibility",[],"Message","");
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",[],"Message","");
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");
output.Message="Optimal Solution Found";
case 1 then
printf("\nMaximum Number of Iterations Exceeded. Output may not be optimal.\n");
output.Message="Maximum Number of Iterations Exceeded. Output may not be optimal";
case 2 then
printf("\nMaximum CPU Time exceeded. Output may not be optimal.\n");
output.Message="Maximum CPU Time exceeded. Output may not be optimal";
case 3 then
printf("\nStop at Tiny Step\n");
output.Message="Stop at Tiny Step";
case 4 then
printf("\nSolved To Acceptable Level\n");
output.Message="Solved To Acceptable Level";
case 5 then
printf("\nConverged to a point of local infeasibility.\n");
output.Message="Converged to a point of local infeasibility";
case 6 then
printf("\nStopping optimization at current point as requested by user.\n");
output.Message="Stopping optimization at current point as requested by user";
case 7 then
printf("\nFeasible point for square problem found.\n");
output.Message="Feasible point for square problem found";
case 8 then
printf("\nIterates diverging; problem might be unbounded.\n");
output.Message="Iterates diverging; problem might be unbounded";
case 9 then
printf("\nRestoration Failed!\n");
output.Message="Restoration Failed!";
case 10 then
printf("\nError in step computation (regularization becomes too large?)!\n");
output.Message="Error in step computation (regularization becomes too large?)!";
case 11 then
printf("\nProblem has too few degrees of freedom.\n");
output.Message="Problem has too few degrees of freedom";
case 12 then
printf("\nInvalid option thrown back by Ipopt\n");
output.Message="Invalid option thrown back by Ipopt";
case 13 then
printf("\nNot enough memory.\n");
output.Message="Not enough memory";
case 15 then
printf("\nINTERNAL ERROR: Unknown SolverReturn value - Notify Ipopt Authors.\n");
output.Message="INTERNAL ERROR: Unknown SolverReturn value - Notify Ipopt Authors";
else
printf("\nInvalid status returned. Notify the Toolbox authors\n");
output.Message="Invalid status returned. Notify the Toolbox authors";
break;
end
endfunction