Go to the documentation of this file.
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.;
Index y_index(Index jx, Index it) const
Translation of mesh point indices to NLP variable indices for y(x_ij)
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 Max(Index a, Index b)
Number beta_
Weighting parameter in PDE.
Class to organize all the data required by the algorithm.
Index Min(Index a, Index b)
bool phi_dydy_always_zero()
Number lb_y_
overall lower bound on y
Class implemented the NLP discretization of.
bool phi_dydy_always_zero()
Class for all IPOPT specific calculated quantities.
bool phi_dydy_always_zero()
Number * a_y_
Array for weighting function a_y in objective.
double Number
Type of all numbers.
virtual bool eval_g(Index n, const Number *x, bool new_x, Index m, Number *g)
Method to return the constraint residuals.
Number ub_u_
overall upper bound on u
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.
Number Number Index Number Number Index Index Index Eval_F_CB Eval_G_CB eval_g
Callback function for evaluating constraint functions.
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...
Number dt_
Step size in t direction.
Number * y_T_
Array for the target profile for y in objective.
Number phi_dydy(Number y)
Base class for parabolic and elliptic control problems, as formulated by Hans Mittelmann as problem (...
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 phi_dydy(Number y)
bool phi_dydy_always_zero()
Number Number Index Number Number Index nele_jac
Number of non-zero elements in constraint Jacobian.
int Index
Type of all indices of vectors, matrices etc.
Number Number Index Number Number Index Index Index Eval_F_CB eval_f
Callback function for evaluating objective function.
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.
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...
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...
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 phi_dydy(Number y)
Index Nt_
Number of mesh points in t direction.
Number ub_y_
overall upper bound on y
Number dx_
Step size in x direction.
Number T_
Upper bound on t.
Number alpha_
Weighting parameter for the control target deviation functional in the objective.
bool phi_dydy_always_zero()
virtual bool eval_f(Index n, const Number *x, bool new_x, Number &obj_value)
Method to return the objective value.
Number phi_dydy(Number y)
Index Nx_
Number of mesh points in x direction.
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 Number * g_scaling
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.
Number t_grid(Index i) const
Compute the grid coordinate for given index in t direction.
Number * a_u_
Array for weighting function a_u in objective.
Number Number * g
Values of constraint at final point (output only - ignored if set to NULL)
IndexStyleEnum
overload this method to return the number of variables and constraints, and the number of non-zeros i...
Number l_
Upper bound on x.
Number x_grid(Index j) const
Compute the grid coordinate for given index in x direction.
Number phi_dydy(Number y)
SolverReturn
enum for the return from the optimize algorithm (obviously we need to add more)
Number lb_u_
overall lower bound on u
MittelmannParaCntrlBase()
Constructor.
virtual ~MittelmannParaCntrlBase()
Default destructor.
virtual bool eval_grad_f(Index n, const Number *x, bool new_x, Number *grad_f)
Method to return the gradient of the objective.
Number Number Index Number Number Index Index Index index_style
indexing style for iRow & jCol, 0 for C style, 1 for Fortran style
Index u_index(Index it) const
virtual bool InitializeProblem(Index N)
Initialize internal parameters, where N is a parameter determining the problme size.
Number Number * x_scaling
Number Number Index m
Number of constraints.