diff options
author | Harpreet | 2015-10-26 15:46:23 +0530 |
---|---|---|
committer | Harpreet | 2015-10-26 15:46:23 +0530 |
commit | 1e82cbeeb97000ce243fb4ef6566d0b554473713 (patch) | |
tree | b9d99d4784c6e6317d67fd9f43baf7f5198a1f53 | |
parent | 4b0ab259b440fa1bb57494dcbdffb865e293b768 (diff) | |
download | FOSSEE-Optimization-toolbox-1e82cbeeb97000ce243fb4ef6566d0b554473713.tar.gz FOSSEE-Optimization-toolbox-1e82cbeeb97000ce243fb4ef6566d0b554473713.tar.bz2 FOSSEE-Optimization-toolbox-1e82cbeeb97000ce243fb4ef6566d0b554473713.zip |
initial guess added
-rw-r--r-- | jar/scilab_en_US_help.jar | bin | 184605 -> 184605 bytes | |||
-rw-r--r-- | macros/lib | bin | 480 -> 480 bytes | |||
-rw-r--r-- | macros/names | 2 | ||||
-rw-r--r-- | macros/qpipopt.bin | bin | 27036 -> 29496 bytes | |||
-rw-r--r-- | macros/qpipopt.sci | 24 | ||||
-rw-r--r-- | macros/qpipoptmat.bin | bin | 0 -> 31280 bytes | |||
-rw-r--r-- | macros/qpipoptmat.sci | 214 | ||||
-rw-r--r-- | macros/qpipoptmat.sci~ (renamed from macros/qpipopt_mat.sci) | 0 | ||||
-rw-r--r-- | sci_gateway/cpp/QuadNLP.hpp | 7 | ||||
-rw-r--r-- | sci_gateway/cpp/QuadNLP.hpp~ | 134 | ||||
-rwxr-xr-x | sci_gateway/cpp/libFAMOS.so | bin | 122834 -> 126936 bytes | |||
-rw-r--r-- | sci_gateway/cpp/sci_QuadNLP.cpp | 4 | ||||
-rw-r--r-- | sci_gateway/cpp/sci_QuadNLP.cpp~ | 253 | ||||
-rw-r--r-- | sci_gateway/cpp/sci_ipopt.cpp | 41 | ||||
-rw-r--r-- | sci_gateway/cpp/sci_ipopt.cpp~ | 409 |
15 files changed, 1074 insertions, 14 deletions
diff --git a/jar/scilab_en_US_help.jar b/jar/scilab_en_US_help.jar Binary files differBinary files differindex 61cfd75..82b1a76 100644 --- a/jar/scilab_en_US_help.jar +++ b/jar/scilab_en_US_help.jar diff --git a/macros/names b/macros/names index ef678fc..40e5934 100644 --- a/macros/names +++ b/macros/names @@ -1,5 +1,5 @@ qpipopt -qpipopt_mat +qpipoptmat setOptions symphony symphony_call diff --git a/macros/qpipopt.bin b/macros/qpipopt.bin Binary files differindex 757b5b6..0cdc0d9 100644 --- a/macros/qpipopt.bin +++ b/macros/qpipopt.bin diff --git a/macros/qpipopt.sci b/macros/qpipopt.sci index 8dfc6c3..a1d094e 100644 --- a/macros/qpipopt.sci +++ b/macros/qpipopt.sci @@ -15,6 +15,7 @@ function [xopt,fopt,exitflag,output,lambda] = qpipopt (varargin) // // Calling Sequence // xopt = qpipopt(nbVar,nbCon,Q,p,LB,UB,conMatrix,conLB,conUB) + // xopt = qpipopt(nbVar,nbCon,Q,p,LB,UB,conMatrix,conLB,conUB,x0) // [xopt,fopt,exitflag,output,lamda] = qpipopt( ... ) // // Parameters @@ -26,7 +27,8 @@ function [xopt,fopt,exitflag,output,lambda] = qpipopt (varargin) // UB : a n x 1 matrix of doubles, where n is number of variables, contains upper bounds of the variables. // conMatrix : a m x n matrix of doubles, where n is number of variables and m is number of constraints, contains matrix representing the constraint matrix // conLB : a m x 1 matrix of doubles, where m is number of constraints, contains lower bounds of the constraints. - // conUB : a m x 1 matrix of doubles, where m is number of constraints, contains upper bounds of the constraints. + // conUB : a m x 1 matrix of doubles, where m is number of constraints, contains upper bounds of the constraints. + // x0 : a m x 1 matrix of doubles, where m is number of constraints, contains initial guess of variables. // xopt : a 1xn matrix of doubles, the computed solution of the optimization problem. // fopt : a 1x1 matrix of doubles, the function value at x. // exitflag : Integer identifying the reason the algorithm terminated. @@ -92,8 +94,8 @@ function [xopt,fopt,exitflag,output,lambda] = qpipopt (varargin) [lhs , rhs] = argn(); //To check the number of argument given by user - if ( rhs ~= 9 ) then - errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 9"), "qpipopt", rhs); + if ( rhs < 9 | rhs > 10 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 9 or 10"), "qpipopt", rhs); error(errmsg) end @@ -108,6 +110,12 @@ function [xopt,fopt,exitflag,output,lambda] = qpipopt (varargin) conLB = varargin(8); conUB = varargin(9); + if ( rhs<10 ) then + x0 = repmat(0,1,nbVar) + else + x0 = varargin(10); + end + //IPOpt wants it in row matrix form p = p'; LB = LB'; @@ -160,11 +168,17 @@ function [xopt,fopt,exitflag,output,lambda] = qpipopt (varargin) //Check the size of constraints of Upper Bound which should equal to the number of constraints if ( size(conUB,2) ~= nbCon) then - errmsg = msprintf(gettext("%s: The Upper Bound of constraints is not equal to the number of constraints"), "qp_ipopt"); + errmsg = msprintf(gettext("%s: The Upper Bound of constraints is not equal to the number of constraints"), "qpipopt"); + error(errmsg); + end + + //Check the size of initial of variables which should equal to the number of variables + if ( size(x0,2) ~= nbVar) then + errmsg = msprintf(gettext("%s: The initial guess of variables is not equal to the number of variables"), "qpipopt"); error(errmsg); end - [xopt,fopt,status,iter,Zl,Zu,lmbda] = solveqp(nbVar,nbCon,Q,p,conMatrix,conLB,conUB,LB,UB); + [xopt,fopt,status,iter,Zl,Zu,lmbda] = solveqp(nbVar,nbCon,Q,p,conMatrix,conLB,conUB,LB,UB,x0); xopt = xopt'; exitflag = status; diff --git a/macros/qpipoptmat.bin b/macros/qpipoptmat.bin Binary files differnew file mode 100644 index 0000000..68c3988 --- /dev/null +++ b/macros/qpipoptmat.bin diff --git a/macros/qpipoptmat.sci b/macros/qpipoptmat.sci new file mode 100644 index 0000000..2f3e911 --- /dev/null +++ b/macros/qpipoptmat.sci @@ -0,0 +1,214 @@ +// Copyright (C) 2015 - IIT Bombay - FOSSEE +// +// Author: Harpreet Singh +// Organization: FOSSEE, IIT Bombay +// Email: harpreet.mertia@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] = qpipoptmat (varargin) + // Solves a linear quadratic problem. + // + // Calling Sequence + // xopt = qpipoptmat(nbVar,nbCon,Q,p,LB,UB,conMatrix,conLB,conUB) + // x = qpipoptmat(H,f) + // x = qpipoptmat(H,f,A,b) + // x = qpipoptmat(H,f,A,b,Aeq,beq) + // x = qpipoptmat(H,f,A,b,Aeq,beq,lb,ub) + // [xopt,fopt,exitflag,output,lamda] = qpipoptmat( ... ) + // + // Parameters + // H : a n x n matrix of doubles, where n is number of variables, represents coefficients of quadratic in the quadratic problem. + // f : a n x 1 matrix of doubles, where n is number of variables, represents coefficients of linear in the quadratic problem + // A : a m x n matrix of doubles, represents the linear coefficients in the inequality constraints + // b : a column vector of doubles, represents the linear coefficients in the inequality constraints + // Aeq : a meq x n matrix of doubles, represents the linear coefficients in the equality constraints + // beq : a vector of doubles, represents the linear coefficients in the equality constraints + // LB : a n x 1 matrix of doubles, where n is number of variables, contains lower bounds of the variables. + // UB : a n x 1 matrix of doubles, where n is number of variables, contains upper bounds of the variables. + // xopt : a nx1 matrix of doubles, the computed solution of the optimization problem. + // fopt : a 1x1 matrix of doubles, the function value at x. + // exitflag : Integer identifying the reason the algorithm terminated. + // output : Structure containing information about the optimization. + // lambda : Structure containing the Lagrange multipliers at the solution x (separated by constraint type). + // + // Description + // Search the minimum of a constrained linear quadratic optimization problem specified by : + // find the minimum of f(x) such that + // + // <latex> + // \begin{eqnarray} + // &\mbox{min}_{x} + // & 1/2*x'*H*x + f'*x \\ + // & \text{subject to} & A.x \leq b \\ + // & & Aeq.x \leq beq \\ + // & & lb \leq x \leq ub \\ + // \end{eqnarray} + // </latex> + // + // We are calling IPOpt for solving the quadratic problem, IPOpt is a library written in C++. The code has been written by Andreas Wächter and Carl Laird. + // + // 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]; + // //and minimize 0.5*x'*Q*x + p'*x with + // f=[1; 2; 3; 4; 5; 6]; H=eye(6,6); + // [xopt,fopt,exitflag,output,lambda]=qpipoptmat(H,f,A,b,Aeq,beq,lb,ub) + // clear H f A b Aeq beq lb ub; + // + // Examples + // //Find the value of x that minimize following function + // // f(x) = 0.5*x1^2 + x2^2 - x1*x2 - 2*x1 - 6*x2 + // // Subject to: + // // x1 + x2 ≤ 2 + // // –x1 + 2x2 ≤ 2 + // // 2x1 + x2 ≤ 3 + // // 0 ≤ x1, 0 ≤ x2. + // 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]; + // [xopt,fopt,exitflag,output,lambda] = qpipoptmat(H,f,A,b,[],[],lb,ub) + // + // Authors + // Keyur Joshi, Saikiran, Iswarya, Harpreet Singh + + +//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 == 3 | rhs == 5 | rhs == 7 | rhs > 8 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be in the set of [2 4 6 8]"), "qpipopt", rhs); + error(errmsg) + end + + H = varargin(1); + f = varargin(2); + nbVar = size(H,1); + + + if ( rhs<2 ) then + A = [] + b = [] + else + A = varargin(3); + b = varargin(4); + end + + if ( rhs<4 ) then + Aeq = [] + beq = [] + else + Aeq = varargin(5); + beq = varargin(6); + end + + if ( rhs<6 ) then + LB = repmat(-%inf,nbVar,1); + UB = repmat(%inf,nbVar,1); + else + LB = varargin(7); + UB = varargin(8); + end + + nbConInEq = size(A,1); + nbConEq = size(Aeq,1); + + //Checking the H matrix which needs to be a symmetric matrix + if ( H~=H') then + errmsg = msprintf(gettext("%s: H is not a symmetric matrix"), "qpipoptmat"); + error(errmsg); + end + + //Check the size of H which should equal to the number of variable + if ( size(H) ~= [nbVar nbVar]) then + errmsg = msprintf(gettext("%s: The Size of H is not equal to the number of variables"), "qpipoptmat"); + 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 Size of f is not equal to the number of variables"), "qpipoptmat"); + error(errmsg); + 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 size of inequality constraints is not equal to the number of variables"), "qpipoptmat"); + 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 size of equality constraints is not equal to the number of variables"), "qpipoptmat"); + 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"), "qpipoptmat"); + 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"), "qpipoptmat"); + 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 Lower Bound of inequality constraints is not equal to the number of constraints"), "qpipoptmat"); + 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 Upper Bound of equality constraints is not equal to the number of constraints"), "qpipoptmat"); + error(errmsg); + end + + //Converting it into ipopt format + f = f'; + LB = LB'; + UB = UB'; + conMatrix = [Aeq;A]; + nbCon = size(conMatrix,1); + conLB = [beq; repmat(-%inf,nbConInEq,1)]'; + conUB = [beq;b]' ; + [xopt,fopt,status,iter,Zl,Zu,lmbda] = solveqp(nbVar,nbCon,H,f,conMatrix,conLB,conUB,LB,UB); + + xopt = xopt'; + exitflag = status; + output = struct("Iterations" , []); + output.Iterations = iter; + lambda = struct("lower" , [], .. + "upper" , [], .. + "ineqlin" , [], .. + "eqlin" , []); + + lambda.lower = Zl; + lambda.upper = Zu; + lambda.eqlin = lmbda(1:nbConEq); + lambda.ineqlin = lmbda(nbConEq+1:nbCon); + + +endfunction diff --git a/macros/qpipopt_mat.sci b/macros/qpipoptmat.sci~ index 4c72216..4c72216 100644 --- a/macros/qpipopt_mat.sci +++ b/macros/qpipoptmat.sci~ diff --git a/sci_gateway/cpp/QuadNLP.hpp b/sci_gateway/cpp/QuadNLP.hpp index 6f01241..ec40195 100644 --- a/sci_gateway/cpp/QuadNLP.hpp +++ b/sci_gateway/cpp/QuadNLP.hpp @@ -47,6 +47,9 @@ class QuadNLP : public TNLP const Number *varLB_= NULL; //varLB_ is a pointer to a matrix of size of 1*numVar_ // with lower bounds of all variables. + + const Number *varGuess_= NULL; //varGuess_ is a pointer to a matrix of size of 1*numVar_ + // with initial guess of all variables. Number *finalX_= NULL; //finalX_ is a pointer to a matrix of size of 1*numVar_ // with final value for the primal variables. @@ -72,8 +75,8 @@ class QuadNLP : public TNLP /* * Constructor */ - QuadNLP(Index nV, Index nC, Number *qM, Number *lM, Number *cM, Number *cUB, Number *cLB, Number *vUB, Number *vLB): - numVars_(nV),numConstr_(nC),qMatrix_(qM),lMatrix_(lM),conMatrix_(cM),conUB_(cUB),conLB_(cLB),varUB_(vUB),varLB_(vLB),finalX_(0), finalZl_(0), finalZu_(0), finalObjVal_(1e20){ } + QuadNLP(Index nV, Index nC, Number *qM, Number *lM, Number *cM, Number *cUB, Number *cLB, Number *vUB, Number *vLB,Number *vG): + numVars_(nV),numConstr_(nC),qMatrix_(qM),lMatrix_(lM),conMatrix_(cM),conUB_(cUB),conLB_(cLB),varUB_(vUB),varLB_(vLB),varGuess_(vG),finalX_(0), finalZl_(0), finalZu_(0), finalObjVal_(1e20){ } /* Go to : diff --git a/sci_gateway/cpp/QuadNLP.hpp~ b/sci_gateway/cpp/QuadNLP.hpp~ new file mode 100644 index 0000000..5e1047f --- /dev/null +++ b/sci_gateway/cpp/QuadNLP.hpp~ @@ -0,0 +1,134 @@ +/* + * Quadratic Programming Toolbox for Scilab using IPOPT library + * Authors : + Sai Kiran + Keyur Joshi + Iswarya + + + * Optimizing (minimizing) the quadratic objective function having any number of variables and linear constraints. + * +*/ + +#ifndef __QuadNLP_HPP__ +#define __QuadNLP_HPP__ + +#include "IpTNLP.hpp" +extern "C"{ +#include <sciprint.h> + +} +using namespace Ipopt; + +class QuadNLP : public TNLP +{ + private: + Index numVars_; // Number of variables. + + Index numConstr_; // Number of constraints. + + const Number *qMatrix_ = NULL; //qMatrix_ is a pointer to matrix of size numVars X numVars_ + // with coefficents of quadratic terms in objective function. + + const Number *lMatrix_ = NULL;//lMatrix_ is a pointer to matrix of size 1*numVars_ + // with coefficents of linear terms in objective function. + + const Number *conMatrix_ = NULL;//conMatrix_ is a pointer to matrix of size numConstr X numVars + // with coefficients of terms in a each objective in each row. + + const Number *conUB_= NULL; //conUB_ is a pointer to a matrix of size of 1*numConstr_ + // with upper bounds of all constraints. + + const Number *conLB_ = NULL; //conLB_ is a pointer to a matrix of size of 1*numConstr_ + // with lower bounds of all constraints. + + const Number *varUB_= NULL; //varUB_ is a pointer to a matrix of size of 1*numVar_ + // with upper bounds of all variables. + + const Number *varLB_= NULL; //varLB_ is a pointer to a matrix of size of 1*numVar_ + // with lower bounds of all variables. + + const Number *varGuess_= NULL; //varGuess_ is a pointer to a matrix of size of 1*numVar_ + // with initial guess of all variables. + + Number *finalX_= NULL; //finalX_ is a pointer to a matrix of size of 1*numVar_ + // with final value for the primal variables. + + Number *finalZl_= NULL; //finalZl_ is a pointer to a matrix of size of 1*numVar_ + // with final values for the lower bound multipliers + + Number *finalZu_= NULL; //finalZu_ is a pointer to a matrix of size of 1*numVar_ + // with final values for the upper bound multipliers + + Number *finalLambda_= NULL; //finalLambda_ is a pointer to a matrix of size of 1*numConstr_ + // with final values for the upper bound multipliers + + Number finalObjVal_; //finalObjVal_ is a scalar with the final value of the objective. + + int iter_; //Number of iteration. + + int status_; //Solver return status + + QuadNLP(const QuadNLP&); + QuadNLP& operator=(const QuadNLP&); + public: + /* + * Constructor + */ + QuadNLP(Index nV, Index nC, Number *qM, Number *lM, Number *cM, Number *cUB, Number *cLB, Number *vUB, Number *vLB): + numVars_(nV),numConstr_(nC),qMatrix_(qM),lMatrix_(lM),conMatrix_(cM),conUB_(cUB),conLB_(cLB),varUB_(vUB),varLB_(vLB),finalX_(0), finalZl_(0), finalZu_(0), finalObjVal_(1e20){ } + + + /* Go to : + + http://www.coin-or.org/Ipopt/documentation/node23.html#SECTION00053130000000000000 + For details about these below methods. + */ + virtual ~QuadNLP(); + virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, + Index& nnz_h_lag, IndexStyleEnum& index_style); + virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u, + Index m, Number* g_l, Number* g_u); + virtual bool get_starting_point(Index n, bool init_x, Number* x, + bool init_z, Number* z_L, Number* z_U, + Index m, bool init_lambda, + Number* lambda); + virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value); + virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f); + virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g); + virtual bool eval_jac_g(Index n, const Number* x, bool new_x, + Index m, Index nele_jac, Index* iRow, Index *jCol, + Number* values); + virtual bool eval_h(Index n, const Number* x, bool new_x, + Number obj_factor, Index m, const Number* lambda, + bool new_lambda, Index nele_hess, Index* iRow, + Index* jCol, Number* values); + virtual void finalize_solution(SolverReturn status, + Index n, const Number* x, const Number* z_L, const Number* z_U, + Index m, const Number* g, const Number* lambda, Number obj_value, + const IpoptData* ip_data, + IpoptCalculatedQuantities* ip_cq); + + const double * getX(); //Returns a pointer to a matrix of size of 1*numVar + // with final value for the primal variables. + + const double * getZu(); //Returns a pointer to a matrix of size of 1*numVars + // with final values for the upper bound multipliers + + const double * getZl(); //Returns a pointer to a matrix of size of 1*numVars + // with final values for the upper bound multipliers + + const double * getLambda(); //Returns a pointer to a matrix of size of 1*numConstr + // with final values for the constraint multipliers + + + double getObjVal(); //Returns the output of the final value of the objective. + + double iterCount(); //Returns the iteration count + + int returnStatus(); //Returns the status count + + +}; + +#endif __QuadNLP_HPP__ diff --git a/sci_gateway/cpp/libFAMOS.so b/sci_gateway/cpp/libFAMOS.so Binary files differindex 5e885f4..d3b1e35 100755 --- a/sci_gateway/cpp/libFAMOS.so +++ b/sci_gateway/cpp/libFAMOS.so diff --git a/sci_gateway/cpp/sci_QuadNLP.cpp b/sci_gateway/cpp/sci_QuadNLP.cpp index ddca7cf..99987a2 100644 --- a/sci_gateway/cpp/sci_QuadNLP.cpp +++ b/sci_gateway/cpp/sci_QuadNLP.cpp @@ -101,8 +101,8 @@ bool QuadNLP::get_starting_point(Index n, bool init_x, Number* x, Index m, bool init_lambda, Number* lambda){ if (init_x == true){ //we need to set initial values for vector x - for (Index var=0;var<n;++var) - x[var]=0.0;//initialize with 0 or we can change. + for (Index var=0;var<n;var++) + x[var]=varGuess_[var];//initialize with 0 or we can change. } if (init_z == true){ //we need to provide initial values for vector bound multipliers diff --git a/sci_gateway/cpp/sci_QuadNLP.cpp~ b/sci_gateway/cpp/sci_QuadNLP.cpp~ new file mode 100644 index 0000000..99987a2 --- /dev/null +++ b/sci_gateway/cpp/sci_QuadNLP.cpp~ @@ -0,0 +1,253 @@ +/* + * Quadratic Programming Toolbox for Scilab using IPOPT library + * Authors : + Sai Kiran + Keyur Joshi + Iswarya + */ + +#include "QuadNLP.hpp" +#include "IpIpoptData.hpp" + +extern "C"{ +#include <api_scilab.h> +#include <Scierror.h> +#include <BOOL.h> +#include <localization.h> +#include <sciprint.h> + + +double x_static,i, *op_obj_x = NULL,*op_obj_value = NULL; + +using namespace Ipopt; + +QuadNLP::~QuadNLP() + { + free(finalX_); + free(finalZl_); + free(finalZu_);} + +//get NLP info such as number of variables,constraints,no.of elements in jacobian and hessian to allocate memory +bool QuadNLP::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, Index& nnz_h_lag, IndexStyleEnum& index_style){ + n=numVars_; // Number of variables + m=numConstr_; // Number of constraints + nnz_jac_g = n*m; // No. of elements in Jacobian of constraints + nnz_h_lag = n*(n+1)/2; // No. of elements in lower traingle of Hessian of the Lagrangian. + index_style=C_STYLE; // Index style of matrices + return true; + } + +//get variable and constraint bound info +bool QuadNLP::get_bounds_info(Index n, Number* x_l, Number* x_u, Index m, Number* g_l, Number* g_u){ + + unsigned int i; + for(i=0;i<n;i++){ + x_l[i]=varLB_[i]; + x_u[i]=varUB_[i]; + } + + for(i=0;i<m;i++){ + g_l[i]=conLB_[i]; + g_u[i]=conUB_[i]; + } + return true; + } + +//get value of objective function at vector x +bool QuadNLP::eval_f(Index n, const Number* x, bool new_x, Number& obj_value){ + unsigned int i,j; + obj_value=0; + + for (i=0;i<=n;i++){ + for (j=0;j<=n;j++){ + obj_value+=0.5*x[i]*x[j]*qMatrix_[n*i+j]; + } + obj_value+=x[i]*lMatrix_[i]; + } + return true; + } + +//get value of gradient of objective function at vector x. +bool QuadNLP::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f){ + unsigned int i,j; + for(i=0;i<n;i++) + { + grad_f[i]=lMatrix_[i]; + for(j=0;j<n;j++) + { + grad_f[i]+=(qMatrix_[n*i+j])*x[j]; + } + } + return true; +} + +//Get the values of constraints at vector x. +bool QuadNLP::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g){ + unsigned int i,j; + for(i=0;i<m;i++) + { + g[i]=0; + for(j=0;j<n;j++) + { + g[i]+=x[j]*conMatrix_[i+j*m]; + } + } + return true; +} + +// This method sets initial values for required vectors . For now we are assuming 0 to all values. +bool QuadNLP::get_starting_point(Index n, bool init_x, Number* x, + bool init_z, Number* z_L, Number* z_U, + Index m, bool init_lambda, + Number* lambda){ + if (init_x == true){ //we need to set initial values for vector x + for (Index var=0;var<n;var++) + x[var]=varGuess_[var];//initialize with 0 or we can change. + } + + if (init_z == true){ //we need to provide initial values for vector bound multipliers + for (Index var=0;var<n;++var){ + z_L[var]=0.0; //initialize with 0 or we can change. + z_U[var]=0.0;//initialize with 0 or we can change. + } + } + + if (init_lambda == true){ //we need to provide initial values for lambda values. + for (Index var=0;var<m;++var){ + lambda[var]=0.0; //initialize with 0 or we can change. + } + } + + return true; + } +/* Return either the sparsity structure of the Jacobian of the constraints, or the values for the Jacobian of the constraints at the point x. + +*/ +bool QuadNLP::eval_jac_g(Index n, const Number* x, bool new_x, + Index m, Index nele_jac, Index* iRow, Index *jCol, + Number* values){ + + //It asked for structure of jacobian. + if (values==NULL){ //Structure of jacobian (full structure) + int index=0; + for (int var=0;var<m;++var)//no. of constraints + for (int flag=0;flag<n;++flag){//no. of variables + iRow[index]=var; + jCol[index]=flag; + index++; + } + } + //It asked for values + else { + int index=0; + for (int var=0;var<m;++var) + for (int flag=0;flag<n;++flag) + values[index++]=conMatrix_[var+flag*m]; + } + return true; + } + +/* + * Return either the sparsity structure of the Hessian of the Lagrangian, + * or the values of the Hessian of the Lagrangian for the given values for + * x,lambda,obj_factor. +*/ +bool QuadNLP::eval_h(Index n, const Number* x, bool new_x, + Number obj_factor, Index m, const Number* lambda, + bool new_lambda, Index nele_hess, Index* iRow, + Index* jCol, Number* values){ + + if (values==NULL){ + Index idx=0; + for (Index row = 0; row < n; row++) { + for (Index col = 0; col <= row; col++) { + iRow[idx] = row; + jCol[idx] = col; + idx++; + } + } + } + else { + Index index=0; + for (Index row=0;row < n;++row){ + for (Index col=0; col <= row; ++col){ + values[index++]=obj_factor*(qMatrix_[n*row+col]); + } + } + } + return true; + } + + +void QuadNLP::finalize_solution(SolverReturn status, + Index n, const Number* x, const Number* z_L, const Number* z_U, + Index m, const Number* g, const Number* lambda, Number obj_value, + const IpoptData* ip_data, + IpoptCalculatedQuantities* ip_cq){ + + finalX_ = (double*)malloc(sizeof(double) * numVars_ * 1); + for (Index i=0; i<n; i++) + { + finalX_[i] = x[i]; + } + + finalZl_ = (double*)malloc(sizeof(double) * numVars_ * 1); + for (Index i=0; i<n; i++) + { + finalZl_[i] = z_L[i]; + } + + finalZu_ = (double*)malloc(sizeof(double) * numVars_ * 1); + for (Index i=0; i<n; i++) + { + finalZu_[i] = z_U[i]; + } + + finalLambda_ = (double*)malloc(sizeof(double) * numConstr_ * 1); + for (Index i=0; i<m; i++) + { + finalLambda_[i] = lambda[i]; + } + + iter_ = ip_data->iter_count(); + finalObjVal_ = obj_value; + status_ = status; + + } + + const double * QuadNLP::getX() + { + return finalX_; + } + + const double * QuadNLP::getZl() + { + return finalZl_; + } + + const double * QuadNLP::getZu() + { + return finalZu_; + } + + const double * QuadNLP::getLambda() + { + return finalLambda_; + } + + double QuadNLP::getObjVal() + { + return finalObjVal_; + } + + double QuadNLP::iterCount() + { + return (double)iter_; + } + + int QuadNLP::returnStatus() + { + return status_; + } + +} diff --git a/sci_gateway/cpp/sci_ipopt.cpp b/sci_gateway/cpp/sci_ipopt.cpp index 06796a9..5837df1 100644 --- a/sci_gateway/cpp/sci_ipopt.cpp +++ b/sci_gateway/cpp/sci_ipopt.cpp @@ -51,13 +51,13 @@ bool readSparse(int arg,int *iRows,int *iCols,int *iNbItem,int** piNbItemRow, in int sci_solveqp(char *fname) { - CheckInputArgument(pvApiCtx, 9, 9); // We need total 9 input arguments. + CheckInputArgument(pvApiCtx, 10, 10); // We need total 10 input arguments. CheckOutputArgument(pvApiCtx, 7, 7); // Error management variable SciErr sciErr; - int retVal=0, *piAddressVarQ = NULL,*piAddressVarP = NULL,*piAddressVarCM = NULL,*piAddressVarCUB = NULL,*piAddressVarCLB = NULL, *piAddressVarLB = NULL,*piAddressVarUB = NULL; - double *QItems=NULL,*PItems=NULL,*ConItems=NULL,*conUB=NULL,*conLB=NULL,*varUB=NULL,*varLB=NULL,x,f,iter; + int retVal=0, *piAddressVarQ = NULL,*piAddressVarP = NULL,*piAddressVarCM = NULL,*piAddressVarCUB = NULL,*piAddressVarCLB = NULL, *piAddressVarLB = NULL,*piAddressVarUB = NULL,*piAddressVarG = NULL; + double *QItems=NULL,*PItems=NULL,*ConItems=NULL,*conUB=NULL,*conLB=NULL,*varUB=NULL,*varLB=NULL,*init_guess = NULL,x,f,iter; static unsigned int nVars = 0,nCons = 0; unsigned int temp1 = 0,temp2 = 0; @@ -260,10 +260,43 @@ int sci_solveqp(char *fname) return 0; } + /* get matrix */ + sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarLB, &temp1,&temp2, &varLB); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + //Initial Value of variables from scilab + /* get Address of inputs */ + sciErr = getVarAddressFromPosition(pvApiCtx, 10, &piAddressVarG); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + /* Check that the first input argument is a real matrix (and not complex) */ + if ( !isDoubleType(pvApiCtx, piAddressVarG) || isVarComplex(pvApiCtx, piAddressVarG) ) + { + Scierror(999, "%s: Wrong type for input argument #%d: A real matrix expected.\n", fname, 10); + return 0; + } + + temp1 = 1; + temp2 = nVars; + /* get matrix */ + sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarG, &temp1,&temp2, &init_guess); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + using namespace Ipopt; - SmartPtr<QuadNLP> Prob = new QuadNLP(nVars,nCons,QItems,PItems,ConItems,conUB,conLB,varUB,varLB); + SmartPtr<QuadNLP> Prob = new QuadNLP(nVars,nCons,QItems,PItems,ConItems,conUB,conLB,varUB,varLB,init_guess); SmartPtr<IpoptApplication> app = IpoptApplicationFactory(); app->RethrowNonIpoptException(true); diff --git a/sci_gateway/cpp/sci_ipopt.cpp~ b/sci_gateway/cpp/sci_ipopt.cpp~ new file mode 100644 index 0000000..8d62b21 --- /dev/null +++ b/sci_gateway/cpp/sci_ipopt.cpp~ @@ -0,0 +1,409 @@ +/* + * Quadratic Programming Toolbox for Scilab using IPOPT library + * Authors : + Sai Kiran + Keyur Joshi + Iswarya + */ + + +#include "sci_iofunc.hpp" +#include "IpIpoptApplication.hpp" +#include "QuadNLP.hpp" + +extern "C"{ +#include <api_scilab.h> +#include <Scierror.h> +#include <BOOL.h> +#include <localization.h> +#include <sciprint.h> + +int j; +double *op_x, *op_obj,*p; + +bool readSparse(int arg,int *iRows,int *iCols,int *iNbItem,int** piNbItemRow, int** piColPos, double** pdblReal){ + SciErr sciErr; + int* piAddr = NULL; + int iType = 0; + int iRet = 0; + sciErr = getVarAddressFromPosition(pvApiCtx, arg, &piAddr); + if(sciErr.iErr) { + printError(&sciErr, 0); + return false; + } + sciprint("\ndone\n"); + if(isSparseType(pvApiCtx, piAddr)){ + sciprint("done\n"); + sciErr =getSparseMatrix(pvApiCtx, piAddr, iRows, iCols, iNbItem, piNbItemRow, piColPos, pdblReal); + if(sciErr.iErr) { + printError(&sciErr, 0); + return false; + } + } + + else { + sciprint("\nSparse matrix required\n"); + return false; + } + return true; + } + +int sci_solveqp(char *fname) +{ + + CheckInputArgument(pvApiCtx, 10, 10); // We need total 10 input arguments. + CheckOutputArgument(pvApiCtx, 7, 7); + + // Error management variable + SciErr sciErr; + int retVal=0, *piAddressVarQ = NULL,*piAddressVarP = NULL,*piAddressVarCM = NULL,*piAddressVarCUB = NULL,*piAddressVarCLB = NULL, *piAddressVarLB = NULL,*piAddressVarUB = NULL,*piAddressVarG = NULL; + double *QItems=NULL,*PItems=NULL,*ConItems=NULL,*conUB=NULL,*conLB=NULL,*varUB=NULL,*varLB=NULL,*init_guess = NULL,x,f,iter; + static unsigned int nVars = 0,nCons = 0; + unsigned int temp1 = 0,temp2 = 0; + + + ////////// Manage the input argument ////////// + + + //Number of Variables + getIntFromScilab(1,&nVars); + + //Number of Constraints + getIntFromScilab(2,&nCons); + + temp1 = nVars; + temp2 = nCons; + + //Q matrix from scilab + /* get Address of inputs */ + sciErr = getVarAddressFromPosition(pvApiCtx, 3, &piAddressVarQ); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + /* Check that the first input argument is a real matrix (and not complex) */ + if ( !isDoubleType(pvApiCtx, piAddressVarQ) || isVarComplex(pvApiCtx, piAddressVarQ) ) + { + Scierror(999, "%s: Wrong type for input argument #%d: A real matrix expected.\n", fname, 3); + return 0; + } + + /* get matrix */ + sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarQ, &temp1, &temp1, &QItems); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + //P matrix from scilab + /* get Address of inputs */ + sciErr = getVarAddressFromPosition(pvApiCtx, 4, &piAddressVarP); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + /* Check that the first input argument is a real matrix (and not complex) */ + if ( !isDoubleType(pvApiCtx, piAddressVarP) || isVarComplex(pvApiCtx, piAddressVarP) ) + { + Scierror(999, "%s: Wrong type for input argument #%d: A real matrix expected.\n", fname, 4); + return 0; + } + + temp1 = 1; + temp2 = nVars; + /* get matrix */ + sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarP, &temp1,&temp2, &PItems); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + if (nCons!=0) + { + //conMatrix matrix from scilab + /* get Address of inputs */ + sciErr = getVarAddressFromPosition(pvApiCtx, 5, &piAddressVarCM); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + /* Check that the first input argument is a real matrix (and not complex) */ + if ( !isDoubleType(pvApiCtx, piAddressVarCM) || isVarComplex(pvApiCtx, piAddressVarCM) ) + { + Scierror(999, "%s: Wrong type for input argument #%d: A real matrix expected.\n", fname, 5); + return 0; + } + temp1 = nCons; + temp2 = nVars; + + /* get matrix */ + sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarCM,&temp1, &temp2, &ConItems); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + + //conLB matrix from scilab + /* get Address of inputs */ + sciErr = getVarAddressFromPosition(pvApiCtx, 6, &piAddressVarCLB); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + /* Check that the first input argument is a real matrix (and not complex) */ + if ( !isDoubleType(pvApiCtx, piAddressVarCLB) || isVarComplex(pvApiCtx, piAddressVarCLB) ) + { + Scierror(999, "%s: Wrong type for input argument #%d: A real matrix expected.\n", fname, 6); + return 0; + } + temp1 = nCons; + temp2 = 1; + + /* get matrix */ + sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarCLB,&temp1, &temp2, &conLB); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + //conUB matrix from scilab + /* get Address of inputs */ + sciErr = getVarAddressFromPosition(pvApiCtx, 7, &piAddressVarCUB); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + /* Check that the first input argument is a real matrix (and not complex) */ + if ( !isDoubleType(pvApiCtx, piAddressVarCUB) || isVarComplex(pvApiCtx, piAddressVarCUB) ) + { + Scierror(999, "%s: Wrong type for input argument #%d: A real matrix expected.\n", fname, 7); + return 0; + } + + temp1 = nCons; + temp2 = 1; + + /* get matrix */ + sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarCUB,&temp1, &temp2, &conUB); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + } + + //varLB matrix from scilab + /* get Address of inputs */ + sciErr = getVarAddressFromPosition(pvApiCtx, 8, &piAddressVarLB); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + /* Check that the first input argument is a real matrix (and not complex) */ + if ( !isDoubleType(pvApiCtx, piAddressVarLB) || isVarComplex(pvApiCtx, piAddressVarLB) ) + { + Scierror(999, "%s: Wrong type for input argument #%d: A real matrix expected.\n", fname, 8); + return 0; + } + temp1 = 1; + temp2 = nVars; + + /* get matrix */ + sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarLB, &temp1,&temp2, &varLB); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + //varUB matrix from scilab + /* get Address of inputs */ + sciErr = getVarAddressFromPosition(pvApiCtx, 9, &piAddressVarUB); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + /* Check that the first input argument is a real matrix (and not complex) */ + if ( !isDoubleType(pvApiCtx, piAddressVarUB) || isVarComplex(pvApiCtx, piAddressVarUB) ) + { + Scierror(999, "%s: Wrong type for input argument #%d: A real matrix expected.\n", fname, 9); + return 0; + } + + temp1 = 1; + temp2 = nVars; + + /* get matrix */ + sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarUB, &temp1,&temp2, &varUB); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + /* get matrix */ + sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarLB, &temp1,&temp2, &varLB); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + //Initial Value of variables from scilab + /* get Address of inputs */ + sciErr = getVarAddressFromPosition(pvApiCtx, 10, &piAddressVarG); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + /* Check that the first input argument is a real matrix (and not complex) */ + if ( !isDoubleType(pvApiCtx, piAddressVarG) || isVarComplex(pvApiCtx, piAddressVarG) ) + { + Scierror(999, "%s: Wrong type for input argument #%d: A real matrix expected.\n", fname, 10); + return 0; + } + + temp1 = 1; + temp2 = nVars; + + /* get matrix */ + sciErr = getMatrixOfDouble(pvApiCtx, piAddressVarG, &temp1,&temp2, &init_guess); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + using namespace Ipopt; + + SmartPtr<QuadNLP> Prob = new QuadNLP(nVars,nCons,QItems,PItems,ConItems,conUB,conLB,varUB,varLB); + SmartPtr<IpoptApplication> app = IpoptApplicationFactory(); + app->RethrowNonIpoptException(true); + + // Change some options + // Note: The following choices are only examples, they might not be + // suitable for your optimization problem. + app->Options()->SetNumericValue("tol", 1e-7); + app->Options()->SetStringValue("mu_strategy", "adaptive"); + + // Indicates whether all equality constraints are linear + app->Options()->SetStringValue("jac_c_constant", "yes"); + // Indicates whether all inequality constraints are linear + app->Options()->SetStringValue("jac_d_constant", "yes"); + // Indicates whether the problem is a quadratic problem + app->Options()->SetStringValue("hessian_constant", "yes"); + + // Initialize the IpoptApplication and process the options + ApplicationReturnStatus status; + status = app->Initialize(); + if (status != Solve_Succeeded) { + sciprint("\n*** Error during initialization!\n"); + return0toScilab(); + return (int) status; + } + // Ask Ipopt to solve the problem + + status = app->OptimizeTNLP(Prob); + + double *fX = Prob->getX(); + double ObjVal = Prob->getObjVal(); + double *Zl = Prob->getZl(); + double *Zu = Prob->getZu(); + double *Lambda = Prob->getLambda(); + double iteration = Prob->iterCount(); + int stats = Prob->returnStatus(); + sciErr = createMatrixOfDouble(pvApiCtx, nbInputArgument(pvApiCtx) + 1, 1, nVars, fX); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + sciErr = createMatrixOfDouble(pvApiCtx, nbInputArgument(pvApiCtx) + 2,1,1,&ObjVal); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + sciErr = createMatrixOfInteger32(pvApiCtx, nbInputArgument(pvApiCtx) + 3,1,1,&stats); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + sciErr = createMatrixOfDouble(pvApiCtx, nbInputArgument(pvApiCtx) + 4,1,1,&iteration); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + sciErr = createMatrixOfDouble(pvApiCtx, nbInputArgument(pvApiCtx) + 5, 1, nVars, Zl); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + sciErr = createMatrixOfDouble(pvApiCtx, nbInputArgument(pvApiCtx) + 6, 1, nVars, Zu); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + sciErr = createMatrixOfDouble(pvApiCtx, nbInputArgument(pvApiCtx) + 7, 1, nCons, Lambda); + if (sciErr.iErr) + { + printError(&sciErr, 0); + return 0; + } + + + AssignOutputVariable(pvApiCtx, 1) = nbInputArgument(pvApiCtx) + 1; + AssignOutputVariable(pvApiCtx, 2) = nbInputArgument(pvApiCtx) + 2; + AssignOutputVariable(pvApiCtx, 3) = nbInputArgument(pvApiCtx) + 3; + AssignOutputVariable(pvApiCtx, 4) = nbInputArgument(pvApiCtx) + 4; + AssignOutputVariable(pvApiCtx, 5) = nbInputArgument(pvApiCtx) + 5; + AssignOutputVariable(pvApiCtx, 6) = nbInputArgument(pvApiCtx) + 6; + AssignOutputVariable(pvApiCtx, 7) = nbInputArgument(pvApiCtx) + 7; + + // As the SmartPtrs go out of scope, the reference count + // will be decremented and the objects will automatically + // be deleted. + + + return 0; + } + +} + +/* +hessian_constan +jacobian _constant + +j_s_d constant : yes +*/ + |