From 89eaf1f5aa8a28122d68dfee0ede409f863f2825 Mon Sep 17 00:00:00 2001 From: jack Date: Fri, 11 Oct 2024 17:50:37 +1000 Subject: [PATCH] Added objective scaling Changed dependencies for cmake sub project support Fixed some memory leaks Fixed printf format strings etc Fixed compatibility with newer compilers --- CMakeLists.txt | 2 +- include/psopt.h | 34 ++---------- src/IPOPT_interface.cxx | 2 +- src/evaluate.cxx | 2 +- src/scaling.cxx | 6 ++- src/setup.cxx | 4 +- src/triplet.cxx | 114 ++++++++++++++++++++-------------------- src/workspace.cxx | 5 +- 8 files changed, 71 insertions(+), 98 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 95eb79bc..b7702a8b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,7 +25,7 @@ PUBLIC $ ) -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) diff --git a/include/psopt.h b/include/psopt.h index bfcb697a..566159f3 100755 --- a/include/psopt.h +++ b/include/psopt.h @@ -30,8 +30,8 @@ e-mail: v.m.becerra@ieee.org using std::numeric_limits; namespace PSOPT { - constexpr double inf = std::numeric_limits::infinity(); - constexpr double pi = 4.0*atan(1.0); + constexpr double inf = std::numeric_limits::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); diff --git a/src/IPOPT_interface.cxx b/src/IPOPT_interface.cxx index 84f1e6db..733296f9 100644 --- a/src/IPOPT_interface.cxx +++ b/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); diff --git a/src/evaluate.cxx b/src/evaluate.cxx index ae71ac43..11062a75 100644 --- a/src/evaluate.cxx +++ b/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; diff --git a/src/scaling.cxx b/src/scaling.cxx index 6f27f97d..25642822 100644 --- a/src/scaling.cxx +++ b/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); } diff --git a/src/setup.cxx b/src/setup.cxx index eb2318f3..8c478c29 100644 --- a/src/setup.cxx +++ b/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; } - diff --git a/src/triplet.cxx b/src/triplet.cxx index 4ec9769d..b216b0c2 100644 --- a/src/triplet.cxx +++ b/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: v.m.becerra@ieee.org - -**********************************************************************************************/ - -#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: v.m.becerra@ieee.org + +**********************************************************************************************/ + +#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; iInsertNonZero( 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; ix0->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; } -