matlabprogram.h
Go to the documentation of this file.
1 // Copyright (C) 2007 Peter Carbonetto. All Rights Reserved.
2 // This code is published under the Eclipse Public License.
3 //
4 // Author: Peter Carbonetto
5 // Dept. of Computer Science
6 // University of British Columbia
7 // May 19, 2007
8 
9 #ifndef INCLUDE_MATLABPROGRAM
10 #define INCLUDE_MATLABPROGRAM
11 
12 #include "array.h"
13 #include "matlabexception.h"
14 #include "matlabfunctionhandle.h"
15 #include "matlabmatrix.h"
16 #include "sparsematrix.h"
17 #include "arrayofmatrices.h"
18 #include "multipliers.h"
19 #include "IpTNLP.hpp"
20 #include "mex.h"
21 
22 using Ipopt::TNLP;
25 using Ipopt::IpoptData;
27 
28 // Class MatlabProgram
29 // -----------------------------------------------------------------
30 class MatlabProgram : public TNLP {
31 public:
32 
33  // The constructor. The input "error" will be set to the appropriate
34  // exception object if the IPOPT solver terminates abnormally.
36  const ArrayOfMatrices& ub, const Matrix& constraintlb,
37  const Matrix& constraintub,
43  const MatlabFunctionHandle& iterFunc, const mxArray* auxData,
47 
48  // The destructor.
49  virtual ~MatlabProgram();
50 
51  // Returns the error generated. If no error was generated, returns a
52  // null pointer.
53  char* geterrormsg() const;
54 
55  // Get the number of number of iterations of IPOPT. Should only be
56  // called after IPOPT has converged to a stationary point.
57  int getnumiterations() const { return numiter; };
58 
59  // Method to return some info about the nonlinear program.
60  virtual bool get_nlp_info (int& numVariables, int& numConstraints,
61  int& sizeOfJ, int& sizeOfH,
62  IndexStyleEnum& indexStyle);
63 
64  // Return the bounds for the problem.
65  virtual bool get_bounds_info (int numVariables, double* lbptr,
66  double* ubptr, int numConstraints,
67  double* clbptr, double* cubptr);
68 
69  // Return the starting point for the algorithm.
70  virtual bool get_starting_point (int numVariables, bool initializeVars,
71  double* variables,
72  bool initializez, double* zl,
73  double* zu, int numConstraints,
74  bool initializeLambda, double* lambda);
75 
76  // Compute the value of the objective.
77  virtual bool eval_f (int numVariables, const double* variables,
78  bool ignoreThis, double& objective);
79 
80  // Compute the gradient of the objective.
81  virtual bool eval_grad_f (int numVariables, const double* variables,
82  bool ignoreThis, double* gradient);
83 
84  // Evaluate the constraint residuals.
85  virtual bool eval_g (int numVariables, const double* variables,
86  bool ignoreThis, int numConstraints,
87  double* constraints);
88 
89  // This method either returns: 1.) The structure of the Jacobian
90  // (if "Jacobian" is zero), or 2.) The values of the Jacobian (if
91  // "Jacobian" is not zero).
92  virtual bool eval_jac_g (int numVariables, const double* variables,
93  bool ignoreThis, int numConstraints,
94  int sizeOfJ, int* rows, int *cols,
95  double* Jacobian);
96 
97  // This method either returns: 1.) the structure of the Hessian of
98  // the Lagrangian (if "Hessian" is zero), or 2.) the values of the
99  // Hessian of the Lagrangian (if "Hesson" is not zero).
100  virtual bool eval_h (int numVariables, const double* variables,
101  bool ignoreThis, double sigma,
102  int numConstraints, const double* multipliers,
103  bool ignoreThisToo, int sizeOfH, int* rows,
104  int* cols, double* Hessian);
105 
106  // This method is called when the algorithm is complete.
107  virtual void finalize_solution (SolverReturn status, int numVariables,
108  const double* variables, const double* zl,
109  const double* zu, int numConstraints,
110  const double* constraints,
111  const double* lambda, double objective,
112  const IpoptData* ip_data,
114 
115  // Intermediate callback method. It is called once per iteration
116  // of the IPOPT algorithm.
117  virtual bool intermediate_callback (AlgorithmMode mode,
118  int iteration, double objective,
119  double inf_pr, double inf_du,
120  double mu, double d_norm,
121  double regularization_ize,
122  double alpha_du, double alpha_pr,
123  int ls_trials,
124  const IpoptData* ip_data,
126 
127 protected:
128  const ArrayOfMatrices& x0; // The initial point.
129  const ArrayOfMatrices& lb; // Lower bounds on the variables.
130  const ArrayOfMatrices& ub; // Upper bounds on the variables.
131  const Matrix& constraintlb; // Lower bounds on the constraints.
132  const Matrix& constraintub; // Upper bounds on the constraints.
133  const mxArray* auxData; // Auxiliary data passed to the
134  // Matlab callback routines.
135  ArrayOfMatrices& xsol; // Storage for the solution.
136  ArrayOfMatrices* x; // Current value of the variables
137  // that's passed to the Matlab
138  // callback routines.
139  Multipliers* initialMultipliers; // The initial point of the
140  // Lagrange multipliers.
141  Multipliers* multipliers; // This is used to store the
142  // value of the Lagrangne
143  // multipliers at the solution.
144  Array<double>* lambda; // Current value of the Lagrange
145  // multipliers that's passed to the
146  // Matlab callback routine for
147  // computing the Hessian.
148  int numiter; // Keeps track of the number of
149  // IPOPT iterations.
150 
151  // Array of inputs to a Matlab callback function.
152  mxArray** prhs;
153  mxArray* lambdarhs;
154 
155  // If true, we don't need to compute the Hessian of the Lagrangian.
157 
158  // These next two members store information about the structure of
159  // the sparse Matlab matrix for the Jacobian of the constraints
160  // and the Hessian of the Lagragian.
163 
164  // The following members specify the Matlab callback routines.
171 
172  // The copy constructor and copy assignment operator are kept
173  // private so that they are not used.
174  MatlabProgram (const MatlabProgram& source);
175  MatlabProgram& operator= (const MatlabProgram& source) { return *this; };
176 
177 private:
178 
179  // Return the value of the objective function at the point "x".
180  double computeObjective (const ArrayOfMatrices& x);
181 
182  // Compute the first derivatives of the objective function at the
183  // point "x" and return their values in "grad".
184  void computeGradient (const ArrayOfMatrices& x, ArrayOfMatrices& grad);
185 
186  // Compute the value of each constraint function at the point "x"
187  // and return the result in second input argument, which is an
188  // array of length equal to the number of constraints.
190 
191  // Compute the derivatives of the Jacobian of the constraints at
192  // the point "x". For this function to work correctly, the Matlab
193  // routine must return the expected sparse structure of the
194  // Jacobian matrix.
195  void computeJacobian (const ArrayOfMatrices& x, double* Jacobian);
196 
197  // Compute the derivatives of the Hessian of the Lagrangian at the
198  // point "x".
199  void computeHessian (const ArrayOfMatrices& x,
200  const Array<double>& lambda,
201  double sigma, double* Hessian);
202 
203  // Call the Matlab routine for computing the Jacobian. The input
204  // arguments to the Matlab routine are the current values of the
205  // optimization variables. The result must be a sparse
206  // matrix. Additionally, there is a second input argument passed
207  // to the Matlab routine. If it is true (1), then we ignore the
208  // actual values entered into the sparse matrix, as long as every
209  // entry that might possibly be a value other than zero at any
210  // time be returned as a non-zero value. It is up to the
211  // programmer to properly deallocate the return value at a later
212  // time by calling mxDestroyArray().
214  bool returnStructureOnly = true);
215 
216  // Call the Matlab routine for computing the Hessian. The input
217  // arguments to the Matlab routine are the current values of the
218  // optimization variables. The result must be a sparse, lower
219  // triangular matrix. Additionally, there is a second input
220  // argument passed to the Matlab routine. If it is true (1), then
221  // we ignore the actual values entered into the sparse matrix, as
222  // long as every entry that might possibly be a value other than
223  // zero at any time be returned as a non-zero value. It is up to
224  // the programmer to properly deallocate the return value at a
225  // later time by calling mxDestroyArray().
226  mxArray* callMatlabHessianRoutine (const ArrayOfMatrices& x,
227  const Array<double>& lambda,
228  bool returnStructureOnly = true,
229  double sigma = 0);
230 };
231 
232 #endif