search

headmyshoulder / odeint-v2

odeint - solving ordinary differential equations in c++ v2
http://headmyshoulder.github.com/odeint-v2/
Other
340 stars 101 forks source link

Added Model Predictive Control (MPC) compatibilities #148

Closed arashium closed 9 years ago

arashium commented 9 years ago

I have added a file integrate_adaptive_intr.hpp in which the adaptive integration can have some forced interruptible time points. This fascility is essential for some decision making processes while control optimizations such as Model Predictive Control. No function is changed but there are some functions cloned and tweaked with new name ending with suffix _interruptable. Calling integrate_adaptive_interruptable is similar to integrate_adaptive except for you should give it an extra function to tell it which points should be forced to the simulation. As an example the function below forces times 0, 0.1, 0.2, 0.3, 0.4 ,0.5, ... to the model simulation:

// called once, and determines the next time to be called while simulation
double interrupt(const state_type &x, const double t)
{
    return t+0.1;
}

The unit test for showing how to use this function is added into example folder. The results shows that these changes work successfully without changing the programming concepts of the current function preserving general case compatiblity. Forced times are bolded below:

*** solve normal *** 0 0 3.88978e-07 7.07049e-08 7.77957e-07 1.48501e-07 2.52836e-06 4.98581e-07 1.04052e-05 2.07393e-06 4.58508e-05 9.16287e-06 0.000205356 4.106e-05 0.000923131 0.000184534 0.00415312 0.000828894 0.0186881 0.0037029 0.0702934 0.0135759 0.121899 0.0229524 0.173504 0.0318573 0.22511 0.0403143 0.278719 0.0486496 0.332329 0.0565498 0.385939 0.0640376 0.44152 0.0713884 0.497101 0.0783417 0.552682 0.0849191 0.61034 0.0913668 0.667998 0.0974532 0.725655 0.103199 0.785508 0.108822 0.84536 0.114119 0.905212 0.119109 0.967391 0.123985 1.02957 0.128568 1.09175 0.132874 1.15641 0.137077 1.22107 0.141017 1.28573 0.14471 1.35304 0.148309 1.42035 0.151674 1.48766 0.15482 1.55781 0.157881 1.62797 0.160734 1.69813 0.163395 1.77135 0.165979 1.84457 0.168381 1.91779 0.170614 1.99433 0.172779 2.07087 0.174785 2.1474 0.176643 2.22755 0.178442 2.30769 0.180102 2.38784 0.181634 2.47191 0.183115 2.55599 0.184477 2.64007 0.185729 2.72846 0.186936 2.81685 0.188041 2.90835 0.189087 2.99985 0.190041 3 0.190043 *** solve interruptable *** 0 0 3.88978e-07 7.07049e-08 7.77957e-07 1.48501e-07 2.52836e-06 4.98581e-07 1.04052e-05 2.07393e-06 4.58508e-05 9.16287e-06 0.000205356 4.106e-05 0.000923131 0.000184534 0.00415312 0.000828894 0.0186881 0.0037029 0.0702934 0.0135759 0.1 0.0190325 0.152388 0.028269 0.2 0.0362538 0.253359 0.0447621 0.3 0.0518364 0.354614 0.0597112 0.4 0.065936 0.455885 0.0732227 0.5 0.0786939 0.55717 0.0854344 0.6 0.0902377 0.65847 0.0964715 0.7 0.100683 0.759787 0.106447 0.8 0.110134 0.861121 0.115462 0.9 0.118686 0.962475 0.123611 1 0.126424 1.06385 0.130975 1.1 0.133426 1.16524 0.137631 1.2 0.139761 1.26666 0.143646 1.3 0.145494 1.3681 0.149082 1.4 0.150681 1.46957 0.153995 1.5 0.155374 1.57106 0.158435 1.6 0.159621 1.67258 0.162448 1.7 0.163463 1.77413 0.166074 1.8 0.16694 1.87571 0.169351 1.9 0.170086 1.97731 0.172312 2 0.172933 2.07895 0.174988 2.1 0.175509 2.18061 0.177406 2.2 0.177839 2.28231 0.17959 2.3 0.179948 2.37959 0.181482 2.4 0.181856 2.48571 0.183347 2.5 0.183583 2.56432 0.184606 2.6 0.185145 2.68887 0.186408 2.7 0.186559 2.75008 0.187215 2.8 0.187838 2.89217 0.188909 2.9 0.188995 2.93522 0.189376 3 0.190043 *** end ***

mariomulansky commented 9 years ago

It seems to me you can already get this behavior using integrate_times. There, you also supply the the time points where the observer should be called, and if you use an adpative stepper, the algorithm automatically makes use of step size control. You won't get observer calls at the intermediate points, but I suppose you don't need them anyways?

arashium commented 9 years ago

Please note that for some applications such as MPC, at each interrupt time a function might be called which makes control decision and influence the process and it is not a merely watching the process.
Does integrate_times interpolate the observer or it add the real points where system is called as well?

mariomulansky commented 9 years ago

when used with a dense_output stepper, integrate_times also does interpolation. But if you use the rk_cash_karp stepper for example, you can be sure that the time points are not interpolated as this stepper cant do interpolation. Anyways, keep in mind that the only way to influence the integration is from the system function. I'm not completely sure what you want to do, so maybe you can provide an example where you actually employ a control so I can figure out how this would be best implemented with odeint?

arashium commented 9 years ago

I had a look at integrate_times. The integrate_times only has constant steps. It cannot run ODE45 with optimised variable steps.

mariomulansky commented 9 years ago

that is not true. there are 3 implementations in https://github.com/headmyshoulder/odeint-v2/blob/master/include/boost/numeric/odeint/integrate/detail/integrate_times.hpp for each case: constant step size, adaptive, and dense output.

arashium commented 9 years ago

Yes. You are right. One thing that is still a question to me is that how can I run system and observer both at normal ODE45 points and force each 0.1*n time? Seems integrate_time forces all manual stop points to the system. But runs the observer only at manual stop points only and not ODE45 chosen points.

void rhs( const double x , double &dxdt , const double t ) { dxdt = 3.0/(2.0_t_t) + x/(2.0*t); cout <<"rhs: "<< t << " " << x << endl;

} void write_cout( const double &x , const double t ) { cout <<"obs: "<< t << " " << x << endl; }

integrate_times( make_controlled( 1E-12 , 1E-12 , stepper_type() ) , rhs , x , times , dt , write_cout );

/* obs: 1 0 rhs: 1 0 rhs: 1.02 0.03 rhs: 1.03 0.0440203 rhs: 1.08 0.113238 rhs: 1.08889 0.124931 rhs: 1.1 0.13954 rhs: 1.1 0.139718 rhs: 1.004 0.006 rhs: 1.006 0.00895977 rhs: 1.016 0.0237156 rhs: 1.01778 0.026315 rhs: 1.02 0.0295565 rhs: 1.02 0.0295583 rhs: 1.00169 0.00252821 rhs: 1.00253 0.00378514 rhs: 1.00674 0.010062 rhs: 1.00749 0.0111736 rhs: 1.00843 0.0125616 rhs: 1.00843 0.0125618 rhs: 1.01068 0.0159027 rhs: 1.01181 0.0175607 rhs: 1.01745 0.0258367 rhs: 1.01845 0.0273007 rhs: 1.0197 0.0291279 rhs: 1.0197 0.0291282 rhs: 1.02196 0.0324141 rhs: 1.02309 0.0340449 rhs: 1.02873 0.0421859 rhs: 1.02973 0.043626 rhs: 1.03098 0.0454235 rhs: 1.03098 0.0454238 rhs: 1.03324 0.0486564 rhs: 1.03437 0.050261 rhs: 1.04 0.0582708 rhs: 1.04101 0.0596879 rhs: 1.04226 0.0614566 rhs: 1.04226 0.0614569 rhs: 1.04462 0.0647927 rhs: 1.04581 0.0664481 rhs: 1.05172 0.0747111 rhs: 1.05277 0.0761727 rhs: 1.05409 0.077997 rhs: 1.05409 0.0779974 rhs: 1.05645 0.0812781 rhs: 1.05763 0.0829062 rhs: 1.06355 0.0910339 rhs: 1.0646 0.0924717 rhs: 1.06591 0.0942663 rhs: 1.06591 0.0942666 rhs: 1.06828 0.0974939 rhs: 1.06946 0.0990958 rhs: 1.07537 0.107092 rhs: 1.07642 0.108507 rhs: 1.07774 0.110273 rhs: 1.07774 0.110273 rhs: 1.08022 0.113607 rhs: 1.08146 0.115261 rhs: 1.08767 0.123519 rhs: 1.08877 0.124979 rhs: 1.09015 0.126802 rhs: 1.09015 0.126803 rhs: 1.09212 0.129403 rhs: 1.09311 0.130695 rhs: 1.09803 0.137148 rhs: 1.09891 0.138291 rhs: 1.1 0.139718 rhs: 1.1 0.139718 obs: 1.1 0.139718 rhs: 1.10256 0.143053 rhs: 1.10384 0.144708 rhs: 1.11024 0.15297 rhs: 1.11138 0.154431 rhs: 1.1128 0.156255 rhs: 1.1128 0.156255 rhs: 1.11536 0.159535 rhs: 1.11664 0.161163 rhs: 1.12303 0.169288 rhs: 1.12417 0.170726 rhs: 1.12559 0.17252 rhs: 1.12559 0.17252 rhs: 1.12815 0.175746 rhs: 1.12943 0.177348 rhs: 1.13583 0.185341 rhs: 1.13697 0.186756 rhs: 1.13839 0.188521 rhs: 1.13839 0.188521 rhs: 1.14108 0.191852 rhs: 1.14242 0.193504 rhs: 1.14913 0.201754 rhs: 1.15033 0.203213 rhs: 1.15182 0.205034 rhs: 1.15182 0.205035 rhs: 1.1545 0.20831 rhs: 1.15585 0.209935 rhs: 1.16256 0.218049 rhs: 1.16375 0.219484 rhs: 1.16524 0.221276 rhs: 1.16524 0.221276 rhs: 1.16793 0.224497 rhs: 1.16927 0.226096 rhs: 1.17599 0.234078 rhs: 1.17718 0.235491 rhs: 1.17867 0.237253 rhs: 1.17867 0.237253 rhs: 1.18149 0.240583 rhs: 1.1829 0.242235 rhs: 1.18995 0.250481 rhs: 1.19121 0.251939 rhs: 1.19277 0.25376 rhs: 1.19277 0.25376 rhs: 1.19422 0.255438 rhs: 1.19494 0.256273 rhs: 1.19855 0.260447 rhs: 1.1992 0.261187 rhs: 1.2 0.262112 rhs: 1.2 0.262112 obs: 1.2 0.262112 rhs: 1.2029 0.265448 rhs: 1.20435 0.267103 rhs: 1.2116 0.275367 rhs: 1.21288 0.276828 rhs: 1.21449 0.278652 rhs: 1.21449 0.278653 rhs: 1.21739 0.281933 rhs: 1.21884 0.283561 rhs: 1.22609 0.291688 rhs: 1.22738 0.293125 rhs: 1.22899 0.29492 rhs: 1.22899 0.29492 rhs: 1.23189 0.298147 rhs: 1.23334 0.299748 rhs: 1.24058 0.307742 rhs: 1.24187 0.309157 rhs: 1.24348 0.310922 rhs: 1.24348 0.310922 rhs: 1.24652 0.314251 rhs: 1.24804 0.315902 rhs: 1.25564 0.324146 rhs: 1.25699 0.325605 rhs: 1.25868 0.327425 rhs: 1.25868 0.327425 rhs: 1.26172 0.330698 rhs: 1.26324 0.332322 rhs: 1.27084 0.34043 rhs: 1.27219 0.341864 rhs: 1.27388 0.343655 rhs: 1.27388 0.343655 rhs: 1.27702 0.346984 rhs: 1.27859 0.348635 rhs: 1.28645 0.35688 rhs: 1.28784 0.358338 rhs: 1.28959 0.360158 rhs: 1.28959 0.360158 rhs: 1.29167 0.362327 rhs: 1.29271 0.363407 rhs: 1.29792 0.368797 rhs: 1.29884 0.369752 rhs: 1.3 0.370945 rhs: 1.3 0.370945 obs: 1.3 0.370945 rhs: 1.30324 0.374279 rhs: 1.30486 0.375934 rhs: 1.31295 0.384193 rhs: 1.31439 0.385654 rhs: 1.31618 0.387478 rhs: 1.31618 0.387478 rhs: 1.31942 0.390757 rhs: 1.32104 0.392384 rhs: 1.32913 0.400507 rhs: 1.33057 0.401944 rhs: 1.33237 0.403738 rhs: 1.33237 0.403738 rhs: 1.3356 0.406963 rhs: 1.33722 0.408564 rhs: 1.34531 0.416556 rhs: 1.34675 0.417969 rhs: 1.34855 0.419734 rhs: 1.34855 0.419734 rhs: 1.35195 0.423066 rhs: 1.35365 0.424719 rhs: 1.36214 0.432971 rhs: 1.36365 0.43443 rhs: 1.36554 0.436252 rhs: 1.36554 0.436252 rhs: 1.36894 0.439528 rhs: 1.37064 0.441154 rhs: 1.37913 0.44927 rhs: 1.38064 0.450706 rhs: 1.38253 0.452498 rhs: 1.38253 0.452498 rhs: 1.38602 0.455812 rhs: 1.38777 0.457457 rhs: 1.39651 0.465665 rhs: 1.39806 0.467117 rhs: 1.4 0.46893 rhs: 1.4 0.46893 obs: 1.4 0.46893 rhs: 1.40351 0.472209 rhs: 1.40527 0.473836 rhs: 1.41406 0.481957 rhs: 1.41562 0.483394 rhs: 1.41757 0.485187 rhs: 1.41757 0.485188 rhs: 1.42121 0.488524 rhs: 1.42303 0.490179 rhs: 1.43212 0.498443 rhs: 1.43373 0.499904 rhs: 1.43575 0.501729 rhs: 1.43575 0.501729 rhs: 1.43939 0.50501 rhs: 1.44121 0.506638 rhs: 1.45029 0.514766 rhs: 1.45191 0.516204 rhs: 1.45393 0.517999 rhs: 1.45393 0.517999 rhs: 1.45769 0.521339 rhs: 1.45957 0.522996 rhs: 1.46898 0.531268 rhs: 1.47065 0.532732 rhs: 1.47274 0.534558 rhs: 1.47274 0.534558 rhs: 1.4765 0.537843 rhs: 1.47838 0.539473 rhs: 1.48779 0.54761 rhs: 1.48946 0.549049 rhs: 1.49155 0.550846 rhs: 1.49155 0.550847 rhs: 1.49324 0.552298 rhs: 1.49408 0.553022 rhs: 1.49831 0.556636 rhs: 1.49906 0.557277 rhs: 1.5 0.558078 rhs: 1.5 0.558078 obs: 1.5 0.558078 rhs: 1.50393 0.561433 rhs: 1.5059 0.563097 rhs: 1.51574 0.571406 rhs: 1.51748 0.572876 rhs: 1.51967 0.574711 rhs: 1.51967 0.574711 rhs: 1.5236 0.57801 rhs: 1.52557 0.579647 rhs: 1.53541 0.587821 rhs: 1.53715 0.589267 rhs: 1.53934 0.591071 rhs: 1.53934 0.591072 rhs: 1.54327 0.594317 rhs: 1.54524 0.595928 rhs: 1.55507 0.603971 rhs: 1.55682 0.605393 rhs: 1.55901 0.607169 rhs: 1.55901 0.60717 rhs: 1.56314 0.610523 rhs: 1.5652 0.612187 rhs: 1.57553 0.620494 rhs: 1.57737 0.621963 rhs: 1.57966 0.623797 rhs: 1.57966 0.623797 rhs: 1.58373 0.627046 rhs: 1.58576 0.628659 rhs: 1.59593 0.636709 rhs: 1.59774 0.638133 rhs: 1.6 0.639911 rhs: 1.6 0.639911 obs: 1.6 0.639911 rhs: 1.60427 0.643271 rhs: 1.60641 0.644938 rhs: 1.6171 0.653261 rhs: 1.619 0.654734 rhs: 1.62137 0.656571 rhs: 1.62137 0.656572 rhs: 1.62565 0.659876 rhs: 1.62779 0.661517 rhs: 1.63847 0.669705 rhs: 1.64037 0.671154 rhs: 1.64275 0.672962 rhs: 1.64275 0.672962 rhs: 1.64718 0.676329 rhs: 1.64939 0.678 rhs: 1.66045 0.686342 rhs: 1.66242 0.687817 rhs: 1.66488 0.689659 rhs: 1.66488 0.689659 rhs: 1.66931 0.692971 rhs: 1.67152 0.694615 rhs: 1.68259 0.702822 rhs: 1.68455 0.704274 rhs: 1.68701 0.706086 rhs: 1.68701 0.706087 rhs: 1.68961 0.707999 rhs: 1.69091 0.708951 rhs: 1.6974 0.713709 rhs: 1.69856 0.714552 rhs: 1.7 0.715605 rhs: 1.7 0.715605 obs: 1.7 0.715605 rhs: 1.70464 0.718991 rhs: 1.70696 0.720671 rhs: 1.71856 0.729059 rhs: 1.72063 0.730542 rhs: 1.7232 0.732394 rhs: 1.7232 0.732395 rhs: 1.72785 0.735725 rhs: 1.73017 0.737379 rhs: 1.74177 0.745631 rhs: 1.74383 0.747092 rhs: 1.74641 0.748914 rhs: 1.74641 0.748914 rhs: 1.75105 0.752192 rhs: 1.75337 0.753819 rhs: 1.76497 0.761942 rhs: 1.76704 0.763379 rhs: 1.76961 0.765173 rhs: 1.76961 0.765173 rhs: 1.7745 0.768567 rhs: 1.77694 0.770251 rhs: 1.78914 0.778658 rhs: 1.79131 0.780145 rhs: 1.79402 0.782002 rhs: 1.79402 0.782002 rhs: 1.79522 0.78282 rhs: 1.79581 0.783229 rhs: 1.7988 0.78527 rhs: 1.79934 0.785632 rhs: 1.8 0.786085 rhs: 1.8 0.786085 obs: 1.8 0.786085 rhs: 1.80501 0.789499 rhs: 1.80752 0.791193 rhs: 1.82004 0.799652 rhs: 1.82227 0.801148 rhs: 1.82505 0.803016 rhs: 1.82505 0.803016 rhs: 1.83006 0.806375 rhs: 1.83257 0.808042 rhs: 1.8451 0.816366 rhs: 1.84732 0.817838 rhs: 1.85011 0.819676 rhs: 1.85011 0.819677 rhs: 1.85512 0.822983 rhs: 1.85762 0.824624 rhs: 1.87015 0.832817 rhs: 1.87238 0.834267 rhs: 1.87516 0.836076 rhs: 1.87516 0.836077 rhs: 1.88013 0.839304 rhs: 1.88261 0.840906 rhs: 1.89503 0.848906 rhs: 1.89724 0.850322 rhs: 1.9 0.852089 rhs: 1.9 0.852089 obs: 1.9 0.852089 rhs: 1.90535 0.855515 rhs: 1.90803 0.857215 rhs: 1.92142 0.865704 rhs: 1.9238 0.867206 rhs: 1.92677 0.86908 rhs: 1.92677 0.869081 rhs: 1.93213 0.872452 rhs: 1.93481 0.874125 rhs: 1.94819 0.882481

...

*/

mariomulansky commented 9 years ago

yes, integrate_times give you only observations at the given points and nothing else. However, I don't see many use cases for your requirements. It seems you are not interested in the observation at those specific time points, but to enforce that the rhs is called at those instances, which somehow is not intended with the integrate routines at all

arashium commented 9 years ago

Some digital control systems need to update the model parameters at scheduled times(for example each 0.2 seconds). They need to manipulate the rhs at the specified time.

headmyshoulder commented 9 years ago

I think you can even call integrate_const with a controlled stepper and put the update of the r.h.s. into the observer. It will call the observer each 0.2 sec for example and run adaptive steps between them. I think you will get very similar results to the output of integrate_adaptive_intr.

On 18.01.2015 20:28, Arash wrote:

Some digital control systems need to update the model parameters at scheduled times(for example each 0.2 seconds). They need to manipulate the rhs at the specified time.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/pull/148#issuecomment-70421698.

arashium commented 9 years ago

Thank you very much. One thing remained as a question for me is that since in real systems, there are boundaries, in x'=Ax+Bu where x is our state space vector, can we put limits for the states?

for example: -2<= x0 <= 5 3<= x1 <=9 -4<= x2 <= 4

x=[x0, x1, x2] T

It is out of potential of rhs and probably must be controlled from the integrate solver.

headmyshoulder commented 9 years ago

On 19.01.2015 04:42, Arash wrote:

Thank you very much. One thing that remained as a question for me is that since in real systems, there are boundaries, in x'=Ax+Bu where x is our state space vector, can we put limits for the states?

for example: -2<= x0 <= 5 3<= x1 <=9 -4<= x2 <= 4

x=[x0, x1, x2] T

It is out of potential of rhs and probably must be controlled from the integrate solver.

This is conceptually not possible with odeint. Of course, you can put the limits into the r.h.s but it will destroy continuity of your ode.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/pull/148#issuecomment-70443110.