diff options
Diffstat (limited to 'macros/intqpipopt.sci')
| -rw-r--r-- | macros/intqpipopt.sci | 461 |
1 files changed, 461 insertions, 0 deletions
diff --git a/macros/intqpipopt.sci b/macros/intqpipopt.sci new file mode 100644 index 0000000..f45d9cb --- /dev/null +++ b/macros/intqpipopt.sci @@ -0,0 +1,461 @@ +// Copyright (C) 2016 - 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: Pranav Deshpande and Akshay Miterani +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in + +function varargout = intqpipopt (varargin) + // Solves a linear quadratic problem. + // + // Calling Sequence + // xopt = intqpipopt(H,f) + // xopt = intqpipopt(H,f,intcon) + // xopt = intqpipopt(H,f,intcon,A,b) + // xopt = intqpipopt(H,f,intcon,A,b,Aeq,beq) + // xopt = intqpipopt(H,f,intcon,A,b,Aeq,beq,lb,ub) + // xopt = intqpipopt(H,f,intcon,A,b,Aeq,beq,lb,ub,x0) + // xopt = intqpipopt(H,f,intcon,A,b,Aeq,beq,lb,ub,x0,options) + // xopt = intqpipopt(H,f,intcon,A,b,Aeq,beq,lb,ub,x0,options,"file_path") + // [xopt,fopt,exitflag,output] = intqpipopt( ... ) + // + // Parameters + // H : a symmetric matrix of double, represents coefficients of quadratic in the quadratic problem. + // f : a vector of double, represents coefficients of linear in the quadratic problem + // intcon : a vector of integers, represents which variables are constrained to be integers + // 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 : a vector of double, contains lower bounds of the variables. + // ub : a vector of double, contains upper bounds of the variables. + // x0 : a vector of double, contains initial guess of variables. + // options : a list containing the parameters to be set. + // file_path : path to bonmin opt file if used. + // xopt : a vector of double, the computed solution of the optimization problem. + // fopt : a double, the value of the function at x. + // exitflag : The exit status. See below for details. + // output : The structure consist of statistics about the optimization. See below for details. + // + // Description + // Search the minimum of a constrained linear quadratic optimization problem specified by : + // + // <latex> + // \begin{eqnarray} + // &\mbox{min}_{x} + // & 1/2⋅x^T⋅H⋅x + f^T⋅x \\ + // & \text{subject to} & A⋅x \leq b \\ + // & & Aeq⋅x = beq \\ + // & & lb \leq x \leq ub \\ + // & & x_i \in \!\, \mathbb{Z}, i \in \!\, intcon\\ + // \end{eqnarray} + // </latex> + // + // The routine calls Bonmin for solving the quadratic problem, Bonmin is a library written in C++. + // + // The exitflag allows to know the status of the optimization which is given back by Bonmin. + // <itemizedlist> + // <listitem>exitflag=0 : Optimal Solution Found. </listitem> + // <listitem>exitflag=1 : InFeasible Solution.</listitem> + // <listitem>exitflag=2 : Output is Continuous Unbounded.</listitem> + // <listitem>exitflag=3 : Limit Exceeded.</listitem> + // <listitem>exitflag=4 : User Interrupt.</listitem> + // <listitem>exitflag=5 : MINLP Error.</listitem> + // </itemizedlist> + // + // For more details on exitflag see the Bonmin page, go to http://www.coin-or.org/Bonmin + // + // The output data structure contains detailed informations about the optimization process. + // It has type "struct" and contains the following fields. + // <itemizedlist> + // <listitem>output.constrviolation: The max-norm of the constraint violation.</listitem> + // </itemizedlist> + // + // Examples + // H = [1 -1; -1 2]; + // f = [-2; -6]; + // A = [1 1; -1 2; 2 1]; + // b = [2; 2; 3]; + // lb=[0,0]; + // ub=[%inf, %inf]; + // intcon = [1 2]; + //[xopt,fopt,status,output]=intqpipopt(H,f,intcon,A,b,[],[],lb,ub) + // + // Examples + // //Find x in R^6 such that: + // Aeq= [1,-1,1,0,3,1; + // -1,0,-3,-4,5,6; + // 2,5,3,0,1,0]; + // beq=[1; 2; 3]; + // A= [0,1,0,1,2,-1; + // -1,0,2,1,1,0]; + // b = [-1; 2.5]; + // lb=[-1000; -10000; 0; -1000; -1000; -1000]; + // ub=[10000; 100; 1.5; 100; 100; 1000]; + // x0 = repmat(0,6,1); + // param = list("MaxIter", 300, "CpuTime", 100); + // //and minimize 0.5*x'*H*x + f'*x with + // f=[1; 2; 3; 4; 5; 6]; H=eye(6,6); + // intcon = [2 4]; + // [xopt,fopt,exitflag,output]=intqpipopt(H,f,intcon,A,b,Aeq,beq,lb,ub,x0,param) + // Authors + // Akshay Miterani and Pranav Deshpande + + //To check the number of input and output argument + [lhs , rhs] = argn(); + + //To check the number of argument given by user + if ( rhs < 2 | rhs == 4 | rhs == 6 | rhs == 8 | rhs > 12 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be in the set of [2 3 5 7 9 10 11 12]"), "intqpipopt", rhs); + error(errmsg) + end + + //To check the number of output arguments + if lhs > 4 then + errmsg = msprintf(gettext("%s: Unexpected number of output arguments: %d provided while should be in the set of [1 2 3 4]"), "intqpipopt", lhs); + end + + H = []; + f = []; + intcon = []; + A = []; + b = []; + Aeq = []; + beq = []; + lb = []; + ub = []; + options = list(); + bonmin_options_file = ''; + + H = varargin(1); + 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"), "intqpipopt"); + error(errmsg); + end + + if ( rhs>=3 ) then + intcon=varargin(3); + end + + if ( rhs<=3 ) then + A = [] + b = [] + else + A = varargin(4); + b = varargin(5); + end + + if ( rhs<6 ) then + Aeq = [] + beq = [] + else + Aeq = varargin(6); + beq = varargin(7); + end + + if ( rhs<8 ) then + lb = repmat(-%inf,nbVar,1); + ub = repmat(%inf,nbVar,1); + else + lb = varargin(8); + ub = varargin(9); + end + + + if ( rhs<10 | size(varargin(10)) ==0 ) then + x0 = repmat(0,nbVar,1) + else + x0 = varargin(10); + end + + if (rhs<11|size(varargin(11))==0) then + options = list(); + else + options = varargin(11); + end + + if(rhs==12) then + bonmin_options_file=varargin(12); + 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 (size(f,2)==0) then + f = repmat(0,nbVar,1); + end + + //Check type of variables + Checktype("intqpipopt", H, "H", 1, "constant") + Checktype("intqpipopt", f, "f", 2, "constant") + Checktype("intqpipopt", intcon, "intcon", 3, "constant") + Checktype("intqpipopt", A, "A", 4, "constant") + Checktype("intqpipopt", b, "b", 5, "constant") + Checktype("intqpipopt", Aeq, "Aeq", 6, "constant") + Checktype("intqpipopt", beq, "beq", 7, "constant") + Checktype("intqpipopt", lb, "lb", 8, "constant") + Checktype("intqpipopt", ub, "ub", 9, "constant") + Checktype("intqpipopt", x0, "x0", 10, "constant") + Checktype("intqpipopt", options, "options", 11, "list") + Checktype("intqpipopt", bonmin_options_file, "bonmin_options_file", 12, "string") + + nbConInEq = size(A,1); + nbConEq = size(Aeq,1); + + // Check if the user gives row vector + // and Changing it to a column matrix + + if (size(f,2)== [nbVar]) then + f=f'; + end + + 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(x0,2)== [nbVar]) then + x0=x0'; + end + + //Checking the H matrix which needs to be a symmetric matrix + if ( ~isequal(H,H')) then + errmsg = msprintf(gettext("%s: H is not a symmetric matrix"), "intqpipopt"); + error(errmsg); + end + + //Check the size of f which should equal to the number of variable + if ( size(f,1) ~= [nbVar]) then + errmsg = msprintf(gettext("%s: The number of rows and columns in H must be equal the number of elements of f"), "intqpipopt"); + error(errmsg); + end + + //Error Checks for intcon + + for i=1:size(intcon,2) + if(intcon(i)>nbVar) then + errmsg = msprintf(gettext("%s: The values inside intcon should be less than the number of variables"), "intqpipopt"); + error(errmsg); + end + + if (intcon(i)<0) then + errmsg = msprintf(gettext("%s: The values inside intcon should be greater than 0 "), "intqpipopt"); + error(errmsg); + end + + if(modulo(intcon(i),1)) then + errmsg = msprintf(gettext("%s: The values inside intcon should be an integer "), "intqpipopt"); + 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 f"), "intqpipopt"); + 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 f"), "intqpipopt"); + 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"), "intqpipopt"); + 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"), "intqpipopt"); + 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,1) ~= 0) then + errmsg = msprintf(gettext("%s: The number of rows in A must be the same as the number of elements of b"), "intqpipopt"); + 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,1) ~= 0) then + errmsg = msprintf(gettext("%s: The number of rows in Aeq must be the same as the number of elements of beq"), "intqpipopt"); + error(errmsg); + end + + //Check the size of initial of variables which should equal to the number of variables + if ( size(x0,1) ~= nbVar) then + warnmsg = msprintf(gettext("%s: Ignoring initial guess of variables as it is not equal to the number of variables"), "intqpipopt"); + warning(warnmsg); + x0 = repmat(0,nbVar,1); + end + + //Check if the user gives a matrix instead of a vector + + if ((size(f,1)~=1)& (size(f,2)~=1)) then + errmsg = msprintf(gettext("%s: f should be a vector"), "intqpipopt"); + error(errmsg); + end + + if (size(lb,1)~=1)& (size(ub,2)~=1) then + errmsg = msprintf(gettext("%s: Lower Bound should be a vector"), "intqpipopt"); + error(errmsg); + end + + if (size(ub,1)~=1)& (size(ub,2)~=1) then + errmsg = msprintf(gettext("%s: Upper Bound should be a vector"), "intqpipopt"); + 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"), "intqpipopt"); + 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"), "intqpipopt"); + error(errmsg); + end + end + + for i = 1:nbConInEq + if (b(i) == -%inf) then + errmsg = msprintf(gettext("%s: Value of b can not be negative infinity"), "intqpipopt"); + error(errmsg); + end + end + + for i = 1:nbConEq + if (beq(i) == -%inf) then + errmsg = msprintf(gettext("%s: Value of beq can not be negative infinity"), "intqpipopt"); + error(errmsg); + end + end + + for i = 1:nbVar + if(lb(i)>ub(i)) then + errmsg = msprintf(gettext("%s: Problem has inconsistent variable bounds"), "intqpipopt"); + error(errmsg); + end + end + + // Checking if the specified options file exists or not + if (~isfile(bonmin_options_file) & ~bonmin_options_file=='') then + error(999, 'The specified options file does not exist!'); + end + + // Validating Options + if (type(options) ~= 15) then + errmsg = msprintf(gettext("%s: Options should be a list "), "cbcintlinprog"); + error(errmsg); + end + + + if (modulo(size(options),2)) then + errmsg = msprintf(gettext("%s: Size of parameters should be even"), "cbcintlinprog"); + error(errmsg); + end + + //Pusing options as required to a double array + optval = []; + ifval =[]; + if length(options) == 0 then + optval = [0 0 0 0 0]; + ifval = [0 0 0 0 0]; + else + optval = [0 0 0 0 0]; + ifval = [0 0 0 0 0]; + for i=1:2:length(options) + select convstr(options(i),'l') + case 'integertolerance' then + ifval(1) = 1; + optval(1) = options(i+1); + case 'maxnodes' then + ifval(2) = 1; + optval(2) = options(i+1); + case 'cputime' then + ifval(3) = 1; + optval(3) = options(i+1); + case 'allowablegap' then + ifval(4) = 1; + optval(4) = options(i+1); + case 'maxiter' then + ifval(5) = 1; + optval(5) = options(i+1); + else + error(999, 'Unknown string argument passed.'); + end + end + end + + //Converting it into bonmin format + f = f'; + lb = lb'; + ub = ub'; + x0 = x0'; + conMatrix = [Aeq;A]; + nbCon = size(conMatrix,1); + conLB = [beq; repmat(-%inf,nbConInEq,1)]'; + conUB = [beq;b]'; + intcon = intcon' + intconSize = length(intcon); + xopt=[] + fopt=[] + status=[] + [xopt,fopt,status] = sci_intqpipopt(nbVar,nbCon,intconSize,H,f,intcon,conMatrix,conLB,conUB,lb,ub,x0,optval,ifval,bonmin_options_file); + xopt = xopt'; + + output.ConstrViolation = max([0;norm(Aeq*xopt-beq, 'inf');(lb'-xopt);(xopt-ub');(A*xopt-b)]); + + varargout(1) = xopt + varargout(2) = fopt + varargout(3) = status + varargout(4) = output + + select status + case 0 then + printf("\nOptimal Solution Found.\n"); + case 1 then + printf("\nInFeasible Solution.\n"); + case 2 then + printf("\nOutput is Continuous Unbounded.s\n"); + case 3 then + printf("\nLimit Exceeded.\n"); + case 4 then + printf("\nUser Interrupt.\n"); + case 5 then + printf("\nMINLP Error.\n"); + else + printf("\nInvalid status returned. Notify the Toolbox authors\n"); + break; + end + +endfunction |