cbcintlinprog added

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