Add objective scaling and other small fixes

CMakeLists.txt

@@ -25,7 +25,7 @@ PUBLIC
     $<INSTALL_INTERFACE:include/${PROJECT_NAME}>
 )
 
-target_link_libraries(${PROJECT_NAME} PRIVATE adolc PkgConfig::ipopt Eigen3::Eigen)
+target_link_libraries(${PROJECT_NAME} PRIVATE adolc Eigen3::Eigen PUBLIC PkgConfig::ipopt)
 add_definitions(-DPSOPT_RELEASE_STRING="${PSOPT_VERSION}" -DHAVE_CSTDDEF)
 
 target_compile_features(${PROJECT_NAME} PUBLIC cxx_delegating_constructors)

include/psopt.h

@@ -30,8 +30,8 @@ e-mail:    [email protected]
 using std::numeric_limits;
 
 namespace PSOPT {
-    constexpr double inf        =   std::numeric_limits<double>::infinity();
-    constexpr double pi         =   4.0*atan(1.0);
+    constexpr double inf    = std::numeric_limits<double>::infinity();
+    const double pi         = 4.0*atan(1.0);
 }
 
 
@@ -192,6 +192,7 @@ struct alg_str {
   string    ipopt_linear_solver;
   double    jac_sparsity_ratio;
   double    hess_sparsity_ratio;
+  double    objective_scaling;
   int       print_level; // 1: detailed output on screen and files (default), 0: no output
   int       save_sparsity_pattern;
   int       nsteps_error_integration;
@@ -1058,35 +1059,6 @@ int get_iphase_offset(Prob& problem, int iphase,Workspace* workspace);
 
 adouble ff_ad(adouble* xad, Workspace* workspace);
 
-
-void rk4_propagate( void (*dae)(adouble* derivatives, adouble* path, adouble* states,
-         adouble* controls, adouble* parameters, adouble& time,
-        adouble* xad, int iphase, Workspace* workspace),
-        MatrixXd& control_trajectory,
-        MatrixXd& time_vector,
-        MatrixXd& initial_state,
-	MatrixXd& parameters,
-        Prob & problem,
-        int iphase,
-        MatrixXd& state_trajectory, Workspace* workspace);
-
-void rkf_propagate( void (*dae)(adouble* derivatives, adouble* path, adouble* states,
-         adouble* controls, adouble* parameters, adouble& time,
-        adouble* xad, int iphase, Workspace* workspace),
-        MatrixXd& control_trajectory,
-        MatrixXd& time_vector,
-        MatrixXd& initial_state,
-	MatrixXd& parameters,
-        double tolerance,
-        double hmin,
-	double hmax,
-        Prob & problem,
-        int iphase,
-        MatrixXd& state_trajectory,
-        MatrixXd& new_time_vector,
-	MatrixXd& new_control_trajectory, Workspace* workspace);
-
-
 void auto_split_observations(Prob& problem, MatrixXd& observation_nodes, MatrixXd& observations);
 
 adouble endpoint_cost_for_parameter_estimation(adouble* initial_states, adouble* final_states, adouble* parameters,adouble& t0, adouble& tf, adouble* xad, int iphase, Workspace* workspace);

src/IPOPT_interface.cxx

@@ -416,7 +416,7 @@ void save_jacobian_sparsity_pattern(Index* rindex, Index* cindex, long nvars, lo
     fprintf(jac_file,"% li %li  %f", ncols, nvars, zero );
 
     for (int i = 1; i< nnz; i++ ) {
-        fprintf(jac_file,"\n%li %li  %f", rindex[i]+1, cindex[i]+1, one ); // SPARSITY PATTERN SAVED USING 1-BASED INDICES (MATLAB-STYLE)
+        fprintf(jac_file,"\n%i %i  %f", rindex[i]+1, cindex[i]+1, one ); // SPARSITY PATTERN SAVED USING 1-BASED INDICES (MATLAB-STYLE)
     }
 
     fclose(jac_file);

src/evaluate.cxx

@@ -259,7 +259,7 @@ void evaluate_solution(Prob& problem,Alg& algorithm,Sol& solution, Workspace* wo
 	psopt_print(workspace,msg);
 	sprintf(msg,"\n_____________________________Statistics per phase______________________________");
 	psopt_print(workspace,msg);
-	sprintf(msg,"\nPhase\t\tNodes\t\tMax ODE Error\tMin ODE error\tMean ODE Error", workspace->current_mesh_refinement_iteration );
+	sprintf(msg,"\nPhase\t\tNodes\t\tMax ODE Error\tMin ODE error\tMean ODE Error");
 	psopt_print(workspace,msg);
 
    solution.mesh_stats[ workspace->current_mesh_refinement_iteration-1 ].epsilon_max = 0;

src/scaling.cxx

@@ -217,8 +217,10 @@ void determine_objective_scaling(MatrixXd& X,Sol& solution, Prob& problem, Alg&
 
   }
 
-
-
+  if (algorithm.objective_scaling != -1) {
+    problem.scale.objective = algorithm.objective_scaling;
+  }
+  printf("Using objective scaling: %e\n", problem.scale.objective);
 }
 
 

src/setup.cxx

@@ -116,7 +116,7 @@ void psopt_level2_setup(Prob& problem, Alg& algorithm)
   problem.bounds.upper.linkage.resize(nlinkages,1);
 
   problem.bounds.lower.linkage.setZero();
-  problem.bounds.lower.linkage.setZero();
+  problem.bounds.upper.linkage.setZero();
 
 
   problem.endpoint_cost               = NULL;
@@ -138,6 +138,7 @@ void psopt_level2_setup(Prob& problem, Alg& algorithm)
   algorithm.nlp_tolerance               = 1.e-6;
   algorithm.jac_sparsity_ratio  	       = 0.5;
   algorithm.hess_sparsity_ratio 	       = 0.2;
+  algorithm.objective_scaling 	        = -1;
   algorithm.hessian                     = "limited-memory";
   algorithm.collocation_method          = "Legendre";
   algorithm.diff_matrix                 = "standard";
@@ -162,4 +163,3 @@ void psopt_level2_setup(Prob& problem, Alg& algorithm)
 
   return;
 }
-

src/triplet.cxx

@@ -1,38 +1,38 @@
-/*********************************************************************************************
-
-This file is part of the PSOPT library, a software tool for computational optimal control
-
-Copyright (C) 2009-2020 Victor M. Becerra
-
-This library is free software; you can redistribute it and/or
-modify it under the terms of the GNU Lesser General Public
-License as published by the Free Software Foundation; either
-version 2.1 of the License, or (at your option) any later version.
-
-This library is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-Lesser General Public License for more details.
-
-You should have received a copy of the GNU Lesser General Public
-License along with this library; if not, write to the Free Software
-Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA,
-or visit http://www.gnu.org/licenses/
-
-Author:    Professor Victor M. Becerra
-Address:   University of Portsmouth
-           School of Energy and Electronic Engineering
-           Portsmouth PO1 3DJ
-           United Kingdom
-e-mail:    [email protected]
-
-**********************************************************************************************/
-
-#include "psopt.h"
-
-// Implementation of TripletSparseMatrix functions
-
-
+/*********************************************************************************************
+
+This file is part of the PSOPT library, a software tool for computational optimal control
+
+Copyright (C) 2009-2020 Victor M. Becerra
+
+This library is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+This library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with this library; if not, write to the Free Software
+Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA,
+or visit http://www.gnu.org/licenses/
+
+Author:    Professor Victor M. Becerra
+Address:   University of Portsmouth
+           School of Energy and Electronic Engineering
+           Portsmouth PO1 3DJ
+           United Kingdom
+e-mail:    [email protected]
+
+**********************************************************************************************/
+
+#include "psopt.h"
+
+// Implementation of TripletSparseMatrix functions
+
+
 TripletSparseMatrix::TripletSparseMatrix(void)
 {
 // Default constructor
@@ -155,11 +155,11 @@ TripletSparseMatrix::TripletSparseMatrix( const TripletSparseMatrix& A) // copy
 TripletSparseMatrix::~TripletSparseMatrix()
 {
     if (a!=NULL)
- 	delete a;
+ 	delete [] a;
     if (RowIndx!= NULL)
-        delete RowIndx;
+        delete [] RowIndx;
     if (ColIndx!= NULL)
-        delete ColIndx;
+        delete [] ColIndx;
 
 }
 
@@ -185,8 +185,8 @@ void TripletSparseMatrix::InsertNonZero(int i, int j, double val)
         }
     }
 
-    if (!eflag) {
-
+    if (!eflag) {
+    	
       anew =  new double[nznew];
       RowNew= new int[nznew];
       ColNew= new int[nznew];
@@ -207,7 +207,7 @@ void TripletSparseMatrix::InsertNonZero(int i, int j, double val)
 
     	RowIndx[nznew-1] = i;
       ColIndx[nznew-1] = j;
-      a[nznew-1]       = val;
+      a[nznew-1]       = val;
 
     }
 
@@ -254,7 +254,7 @@ TripletSparseMatrix& TripletSparseMatrix::operator -= (const TripletSparseMatrix
 
     for (i=0; i<rval.nz; i++)
     {
-        if ( (*this)(rval.RowIndx[i],rval.ColIndx[i])!=0.0 )
+        if ( (*this)(rval.RowIndx[i],rval.ColIndx[i])!=0.0 )
              (*this)(rval.RowIndx[i],rval.ColIndx[i]) -= rval.a[i];
         else
                this->InsertNonZero( rval.RowIndx[i], rval.ColIndx[i], -rval.a[i] );
@@ -301,8 +301,8 @@ TripletSparseMatrix operator *(double Arg, const TripletSparseMatrix& A)
 {
    return (A*Arg);
 }
-
-
+
+
 TripletSparseMatrix TripletSparseMatrix::operator* (const MatrixXd& A) const
 {
   int k;
@@ -329,7 +329,7 @@ TripletSparseMatrix TripletSparseMatrix::operator* (const MatrixXd& A) const
 
   return (sp);
 
-}
+}
 
 
 TripletSparseMatrix elemProduct(const TripletSparseMatrix A, const TripletSparseMatrix& B)
@@ -402,12 +402,12 @@ double& TripletSparseMatrix::operator() (int i,int j)
 
     int ii;
 
-    int eflag = 0;
-
-    if (i>=n || i<0 || j >=m || j<0 ) {
-
-         sp_error_message("Index out of range in operator TripletSparseMatrix::operator(int, int)");    
-
+    int eflag = 0;
+    
+    if (i>=n || i<0 || j >=m || j<0 ) {
+
+         sp_error_message("Index out of range in operator TripletSparseMatrix::operator(int, int)");    
+    
     } 
 
     for (k=0; k< nz; k++)
@@ -471,15 +471,15 @@ void TripletSparseMatrix::Print(const char* text)
    int i;
 
    fprintf(stderr,"\nSparse matrix %s",text);
-   fprintf(stderr,"\nNumber of rows: %li", n);
-   fprintf(stderr,"\nNumber of columns: %li", m);
-   fprintf(stderr,"\nNumber of non-zero elements: %li", nz);
+   fprintf(stderr,"\nNumber of rows: %d", n);
+   fprintf(stderr,"\nNumber of columns: %d", m);
+   fprintf(stderr,"\nNumber of non-zero elements: %d", nz);
    if (n*m!=0) fprintf(stderr,"\nDensity: %lf%%", (  ((double) nz*100)/(double) (n*m)  ) );
    if (nz>0) fprintf(stderr,"\n(Row,Col)\tValue");
 
    for (i=0; i<nz; i++)
    {
-       fprintf(stderr,"\n(%li,%li)\t\t%e", RowIndx[i], ColIndx[i], a[i]);
+       fprintf(stderr,"\n(%d,%d)\t\t%e", RowIndx[i], ColIndx[i], a[i]);
 
    }
    fprintf(stderr,"\n");
@@ -619,7 +619,7 @@ void TripletSparseMatrix::Save(const char* fname) const
 //   fprintf(fp,"%li\t%li\t%li\n", n, m, nz);
 
    for (k=0;k<nz;k++) {
-            fprintf(fp,"%li\t%li\t%e\n",  RowIndx[k]+1, ColIndx[k]+1, a[k] );
+            fprintf(fp,"%d\t%d\t%e\n",  RowIndx[k]+1, ColIndx[k]+1, a[k] );
    }
 
    fclose(fp);
@@ -705,4 +705,4 @@ void TripletSparseMatrix::Transpose()
 void sp_error_message(const char *error_text)
 {
      error_message( error_text );
-}
+}

src/workspace.cxx

@@ -337,8 +337,8 @@ void resize_workspace_vars(Prob& problem, Alg& algorithm, Sol& solution, Workspa
   workspace->x0->resize(nvars,1);
   workspace->lambda->resize(nlp_ncons,1);
 
-  workspace->xlb->resize(nvars,1);
-  workspace->xub->resize(nvars,1);
+  *workspace->xlb = MatrixXd::Zero(nvars,1);
+  *workspace->xub = MatrixXd::Zero(nvars,1);
 
   workspace->nphases = problem.nphases;
 
@@ -535,4 +535,3 @@ work_str::~work_str()
   delete this->grw;
 
 }
-