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,
60 Index m,
bool init_lambda,
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()) {
262 index_style = C_STYLE;
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.;
355 obj_scaling = 1./
Min(dx_,dt_);
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.;