9 #ifndef __MITTELMANNPARACNTRL_HPP__
10 #define __MITTELMANNPARACNTRL_HPP__
17 #include "configall_system.h"
26 # error "don't have header file for math"
30 using namespace Ipopt;
58 virtual bool get_starting_point(
Index n,
bool init_x,
Number*
x,
93 bool& use_x_scaling,
Index n,
95 bool& use_g_scaling,
Index m,
110 virtual bool InitializeProblem(
Index N);
169 return jx + (Nx_+1)*it;
173 return (Nt_+1)*(Nx_+1) + it - 1;
208 typename T::ProblemSpecs p;
211 printf(
"N has to be at least 1.");
229 for (
Index j=0; j<=Nx_; j++) {
230 y_T_[j] = p.y_T(x_grid(j));
233 for (
Index i=1; i<=Nt_; i++) {
234 a_y_[i-1] = p.a_y(t_grid(i));
237 for (
Index i=1; i<=Nt_; i++) {
238 a_u_[i-1] = p.a_u(t_grid(i));
249 typename T::ProblemSpecs p;
251 n = (Nt_+1)*(Nx_+1) + Nt_;
253 m = Nt_*(Nx_-1) + Nt_ + Nt_;
255 nnz_jac_g = 6*Nt_*(Nx_-1) + 3*Nt_ + 4*Nt_;
257 nnz_h_lag = Nx_+1 + Nt_;
258 if (!p.phi_dydy_always_zero()) {
272 typename T::ProblemSpecs p;
275 for (
Index jx=0; jx<=Nx_; jx++) {
276 for (
Index it=1; it<=Nt_; it++) {
277 Index iy = y_index(jx,it);
282 for (
Index i=1; i<=Nt_; i++) {
283 Index iu = u_index(i);
296 for (
Index jx=0; jx<=Nx_; jx++) {
297 Index iy = y_index(jx,0);
298 x_u[iy] = x_l[iy] = p.a(x_grid(jx));
302 for (
Index i=0; i<Nt_*(Nx_-1) + Nt_; i++) {
308 for (
Index i=0; i<Nt_; i++) {
309 g_l[Nt_*(Nx_-1) + Nt_ + i]
310 = g_u[Nt_*(Nx_-1) + Nt_ + i]
321 Index m,
bool init_lambda,
324 DBG_ASSERT(init_x==
true && init_z==
false && init_lambda==
false);
327 for (
Index jx=0; jx<=Nx_; jx++) {
328 for (
Index it=0; it<=Nt_; it++) {
329 x[y_index(jx,it)] = 0.;
333 for (
Index i=1; i<=Nt_; i++) {
334 x[u_index(i)] = (ub_u_+lb_u_)/2.;
356 use_x_scaling =
false;
357 use_g_scaling =
false;
364 bool new_x,
Number& obj_value)
367 Number a =
x[y_index(0,Nt_)] - y_T_[0];
369 for (
Index jx=1; jx<Nx_; jx++) {
370 a =
x[y_index(jx,Nt_)] - y_T_[jx];
373 a =
x[y_index(Nx_,Nt_)] - y_T_[Nx_];
375 obj_value = 0.5*dx_*sum;
379 sum = 0.5*
x[u_index(Nt_)]*
x[u_index(Nt_)];
380 for (
Index it=1; it < Nt_; it++) {
381 sum +=
x[u_index(it)]*
x[u_index(it)];
383 obj_value += 0.5*alpha_*dt_*sum;
388 for (
Index it=1; it<Nt_; it++) {
389 sum += a_y_[it-1]*
x[y_index(Nx_,it)] + a_u_[it-1]*
x[u_index(it)];
391 sum += 0.5*(a_y_[Nt_-1]*
x[y_index(Nx_,Nt_)] + a_u_[Nt_-1]*
x[u_index(Nt_)]);
392 obj_value += dt_*sum;
402 for (
Index jx=0; jx<=Nx_; jx++) {
403 for (
Index it=0; it<=Nt_; it++) {
404 grad_f[y_index(jx,it)] = 0.;
409 grad_f[y_index(0,Nt_)] = 0.5*dx_*(
x[y_index(0,Nt_)] - y_T_[0]);
410 for (
Index jx=1; jx<Nx_; jx++) {
411 grad_f[y_index(jx,Nt_)] = dx_*(
x[y_index(jx,Nt_)] - y_T_[jx]);
413 grad_f[y_index(Nx_,Nt_)] = 0.5*dx_*(
x[y_index(Nx_,Nt_)] - y_T_[Nx_]);
416 for (
Index it=1; it<Nt_; it++) {
417 grad_f[y_index(Nx_,it)] = dt_*a_y_[it-1];
419 grad_f[y_index(Nx_,Nt_)] += 0.5*dt_*a_y_[Nt_-1];
422 for (
Index it=1; it<Nt_; it++) {
423 grad_f[u_index(it)] = alpha_*dt_*
x[u_index(it)] + dt_*a_u_[it-1];
425 grad_f[u_index(Nt_)] = 0.5*dt_*(alpha_*
x[u_index(Nt_)] + a_u_[Nt_-1]);
434 typename T::ProblemSpecs p;
438 Number f = 1./(2.*dx_*dx_);
439 for (
Index jx=1; jx<Nx_; jx++) {
440 for (
Index it=0; it<Nt_; it++) {
441 g[ig] = (
x[y_index(jx,it)]-
x[y_index(jx,it+1)])/dt_
442 + f*(
x[y_index(jx-1,it)] - 2.*
x[y_index(jx,it)]
443 +
x[y_index(jx+1,it)] +
x[y_index(jx-1,it+1)]
444 - 2.*
x[y_index(jx,it+1)] +
x[y_index(jx+1,it+1)]);
449 for (
Index it=1; it<=Nt_; it++) {
450 g[ig] = (
x[y_index(2,it)] - 4.*
x[y_index(1,it)] + 3.*
x[y_index(0,it)]);
455 for (
Index it=1; it<=Nt_; it++) {
457 f*(
x[y_index(Nx_-2,it)] - 4.*
x[y_index(Nx_-1,it)]
458 + 3.*
x[y_index(Nx_,it)]) + beta_*
x[y_index(Nx_,it)]
459 -
x[u_index(it)] + p.phi(
x[y_index(Nx_,it)]);
474 typename T::ProblemSpecs p;
478 if (values == NULL) {
480 for (
Index jx=1; jx<Nx_; jx++) {
481 for (
Index it=0; it<Nt_; it++) {
483 jCol[ijac] = y_index(jx-1,it);
486 jCol[ijac] = y_index(jx,it);
489 jCol[ijac] = y_index(jx+1,it);
492 jCol[ijac] = y_index(jx-1,it+1);
495 jCol[ijac] = y_index(jx,it+1);
498 jCol[ijac] = y_index(jx+1,it+1);
505 for (
Index it=1; it<=Nt_; it++) {
507 jCol[ijac] = y_index(0,it);
510 jCol[ijac] = y_index(1,it);
513 jCol[ijac] = y_index(2,it);
519 for (
Index it=1; it<=Nt_; it++) {
521 jCol[ijac] = y_index(Nx_-2,it);
524 jCol[ijac] = y_index(Nx_-1,it);
527 jCol[ijac] = y_index(Nx_,it);
530 jCol[ijac] = u_index(it);
538 Number f = 1./(2.*dx_*dx_);
540 for (
Index jx=1; jx<Nx_; jx++) {
541 for (
Index it=0; it<Nt_; it++) {
543 values[ijac++] = f2 - 2.*f;
546 values[ijac++] = -f2 - 2.*f;
551 for (
Index it=1; it<=Nt_; it++) {
553 values[ijac++] = -4.;
558 for (
Index it=1; it<=Nt_; it++) {
560 values[ijac++] = -4.*f;
561 values[ijac++] = 3.*f + beta_ + p.phi_dy(
x[y_index(Nx_,it)]);
562 values[ijac++] = -1.;
577 typename T::ProblemSpecs p;
581 if (values == NULL) {
583 for (
Index jx=0; jx<= Nx_; jx++) {
584 iRow[ihes] = y_index(jx,Nt_);
585 jCol[ihes] = y_index(jx,Nt_);
589 for (
Index it=1; it<=Nt_; it++) {
590 iRow[ihes] = u_index(it);
591 jCol[ihes] = u_index(it);
596 if (!p.phi_dydy_always_zero()) {
597 for (
Index it=1; it<=Nt_; it++) {
598 iRow[ihes] = y_index(Nx_,it);
599 jCol[ihes] = y_index(Nx_,it);
606 values[ihes++] = obj_factor*0.5*dx_;
607 for (
Index jx=1; jx<Nx_; jx++) {
608 values[ihes++] = obj_factor*dx_;
610 values[ihes++] = obj_factor*0.5*dx_;
612 for (
Index it=1; it<Nt_; it++) {
613 values[ihes++] = obj_factor*alpha_*dt_;
615 values[ihes++] = obj_factor*0.5*alpha_*dt_;
618 if (!p.phi_dydy_always_zero()) {
619 Index ig = (Nx_-1)*Nt_ + Nt_;
620 for (
Index it=1; it<=Nt_; it++) {
621 values[ihes++] = lambda[ig++]*p.phi_dydy(
x[y_index(Nx_,it)]);
712 return y*pow(fabs(y),3);
716 return 4.*pow(fabs(y),3);
1039 return exp(-4.*t)/4.
1044 return -y*sin(y/10.);
1048 return -y*cos(y/10.)/10. - sin(y/10.);
1052 return y*sin(y/10.)/100.;
Number * x
Input: Starting point Output: Optimal solution.
Number Number Index Number Number Index Index nele_hess
Number of non-zero elements in Hessian of Lagrangian.
Number Number Index Number Number Index Index Index Eval_F_CB Eval_G_CB Eval_Grad_F_CB Eval_Jac_G_CB eval_jac_g
Callback function for evaluating Jacobian of constraint functions.
Number Number * g
Values of constraint at final point (output only - ignored if set to NULL)
Number Number Index Number Number Index nele_jac
Number of non-zero elements in constraint Jacobian.
Number Number Index Number Number Index Index Index Eval_F_CB eval_f
Callback function for evaluating objective function.
Number Number Index Number Number Index Index Index Eval_F_CB Eval_G_CB eval_g
Callback function for evaluating constraint functions.
Number Number Index Number Number Index Index Index Eval_F_CB Eval_G_CB Eval_Grad_F_CB Eval_Jac_G_CB Eval_H_CB eval_h
Callback function for evaluating Hessian of Lagrangian function.
Number Number * x_scaling
Number Number Number * g_scaling
Number Number Index Number Number Index Index Index Eval_F_CB Eval_G_CB Eval_Grad_F_CB eval_grad_f
Callback function for evaluating gradient of objective function.
Number Number Index m
Number of constraints.
Number Number Index Number Number Index Index Index index_style
indexing style for iRow & jCol, 0 for C style, 1 for Fortran style
Class for all IPOPT specific calculated quantities.
Class to organize all the data required by the algorithm.
IndexStyleEnum
overload this method to return the number of variables and constraints, and the number of non-zeros i...
Number phi_dydy(Number y)
bool phi_dydy_always_zero()
Number phi_dydy(Number y)
bool phi_dydy_always_zero()
bool phi_dydy_always_zero()
Number phi_dydy(Number y)
Number phi_dydy(Number y)
bool phi_dydy_always_zero()
Number phi_dydy(Number y)
bool phi_dydy_always_zero()
Base class for parabolic and elliptic control problems, as formulated by Hans Mittelmann as problem (...
Number ub_u_
overall upper bound on u
virtual bool eval_jac_g(Index n, const Number *x, bool new_x, Index m, Index nele_jac, Index *iRow, Index *jCol, Number *values)
Method to return: 1) The structure of the jacobian (if "values" is NULL) 2) The values of the jacobia...
virtual ~MittelmannParaCntrlBase()
Default destructor.
Number alpha_
Weighting parameter for the control target deviation functional in the objective.
virtual bool get_bounds_info(Index n, Number *x_l, Number *x_u, Index m, Number *g_l, Number *g_u)
Method to return the bounds for my problem.
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)
This method is called after the optimization, and could write an output file with the optimal profile...
Number t_grid(Index i) const
Compute the grid coordinate for given index in t direction.
virtual bool eval_grad_f(Index n, const Number *x, bool new_x, Number *grad_f)
Method to return the gradient of the objective.
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)
Method to return the starting point for the algorithm.
MittelmannParaCntrlBase & operator=(const MittelmannParaCntrlBase< T > &)
virtual bool get_scaling_parameters(Number &obj_scaling, bool &use_x_scaling, Index n, Number *x_scaling, bool &use_g_scaling, Index m, Number *g_scaling)
Method for returning scaling parameters.
Number x_grid(Index j) const
Compute the grid coordinate for given index in x direction.
Number ub_y_
overall upper bound on y
Number T_
Upper bound on t.
Number lb_y_
overall lower bound on y
Number beta_
Weighting parameter in PDE.
Index u_index(Index it) const
MittelmannParaCntrlBase()
Constructor.
MittelmannParaCntrlBase(const MittelmannParaCntrlBase< T > &)
Number l_
Upper bound on x.
virtual bool get_nlp_info(Index &n, Index &m, Index &nnz_jac_g, Index &nnz_h_lag, IndexStyleEnum &index_style)
Method to return some info about the nlp.
Index Nt_
Number of mesh points in t direction.
Index Nx_
Number of mesh points in x direction.
virtual bool eval_f(Index n, const Number *x, bool new_x, Number &obj_value)
Method to return the objective value.
Number * a_u_
Array for weighting function a_u in objective.
Number lb_u_
overall lower bound on u
virtual bool eval_g(Index n, const Number *x, bool new_x, Index m, Number *g)
Method to return the constraint residuals.
Index y_index(Index jx, Index it) const
Translation of mesh point indices to NLP variable indices for y(x_ij)
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)
Method to return: 1) The structure of the hessian of the lagrangian (if "values" is NULL) 2) The valu...
virtual bool InitializeProblem(Index N)
Initialize internal parameters, where N is a parameter determining the problme size.
Number * y_T_
Array for the target profile for y in objective.
Number dx_
Step size in x direction.
Number dt_
Step size in t direction.
Number * a_y_
Array for weighting function a_y in objective.
Class implemented the NLP discretization of.
SolverReturn
enum for the return from the optimize algorithm (obviously we need to add more)
Index Max(Index a, Index b)
int Index
Type of all indices of vectors, matrices etc.
Index Min(Index a, Index b)
double Number
Type of all numbers.