31 #ifndef _adevs_hybrid_h_ 32 #define _adevs_hybrid_h_ 35 #include "adevs_models.h" 50 N(N_vars),M(M_event_funcs){}
56 virtual void init(
double* q) = 0;
58 virtual void der_func(
const double* q,
double* dq) = 0;
77 const bool* state_event) = 0;
85 virtual void output_func(
const double *q,
const bool* state_event,
120 double err_tol = 1E-10,
int max_iters = 30,
double alpha = -1.0):
123 max_iters(max_iters),
130 atmp =
new double[A];
132 f[1] =
new double[A];
133 f[0] =
new double[A];
143 virtual void init(
double* q,
double* a) = 0;
150 virtual void alg_func(
const double* q,
const double* a,
double* af) = 0;
155 virtual void der_func(
const double* q,
const double* a,
double* dq) = 0;
159 virtual void state_event_func(
const double* q,
const double* a,
double* z) = 0;
167 virtual void postStep(
double* q,
double* a) = 0;
174 const bool* state_event) = 0;
177 double e,
const Bag<X>& xb) = 0;
180 const bool* state_event,
const Bag<X>& xb) = 0;
182 virtual void output_func(
const double *q,
const double* a,
183 const bool* state_event,
Bag<X>& yb) = 0;
274 void solve(
const double* q);
276 const int A, max_iters;
277 const double err_tol, alpha;
290 template <
typename X>
293 int iter_count = 0, alt, good;
294 double prev_err, err = 0.0, ee, beta, g2, alpha_tmp = alpha;
300 _adevs_dae_se_1_system_solve_try_it_again:
305 alg_func(q,a,f[alt]);
306 for (
int i = 0; i < A; i++)
314 a[i] += alpha_tmp*d[i];
318 while (iter_count < max_iters)
323 alg_func(q,a,f[good]);
325 for (
int i = 0; i < A; i++)
330 ee = fabs(f[good][i]);
331 if (ee > err) err = ee;
334 if (err < err_tol)
return;
339 for (
int i = 0; i < A; i++)
342 if (alpha_tmp < 0.0) alpha_tmp = -alpha_tmp;
343 else alpha_tmp *= -0.5;
344 goto _adevs_dae_se_1_system_solve_try_it_again;
350 for (
int i = 0; i < A; i++)
351 g2 += f[alt][i]*f[alt][i];
352 for (
int i = 0; i < A; i++)
353 beta += f[good][i]*(f[good][i]-f[alt][i]);
356 for (
int i = 0; i < A; i++)
358 d[i] = beta*d[i]-f[good][i];
360 a[i] += alpha_tmp*d[i];
390 virtual double integrate(
double* q,
double h_lim) = 0;
394 virtual void advance(
double* q,
double h) = 0;
424 virtual bool find_events(
bool* events,
const double *qstart,
440 template <
typename X,
class T =
double>
class Hybrid:
450 sys(sys),solver(solver),event_finder(event_finder),
453 q =
new double[sys->
numVars()];
454 q_trial =
new double[sys->
numVars()];
456 event_exists =
false;
458 for (
int i = 0; i < sys->
numVars(); i++) q[i] = q_trial[i];
475 if (!missedOutput.empty())
477 missedOutput.clear();
482 event_happened = event_exists;
485 sys->internal_event(q_trial,event);
489 for (
int i = 0; i < sys->numVars(); i++) q[i] = q_trial[i];
498 bool state_event_exists =
false;
499 event_happened =
true;
503 for (
int i = 0; i < sys->numVars(); i++)
505 solver->advance(q_trial,e);
507 event_finder->find_events(event,q,q_trial,solver,e);
509 if (state_event_exists)
512 sys->confluent_event(q_trial,event,xb);
513 for (
int i = 0; i < sys->numVars(); i++)
517 if (!state_event_exists)
519 solver->advance(q,e);
523 sys->external_event(q,e+e_accum,xb);
527 for (
int i = 0; i < sys->numVars(); i++) q_trial[i] = q[i];
536 if (!missedOutput.empty())
538 missedOutput.clear();
539 if (sigma > 0.0) event_exists =
false;
542 event_happened =
true;
544 sys->confluent_event(q_trial,event,xb);
545 else sys->external_event(q_trial,e_accum+ta(),xb);
548 for (
int i = 0; i < sys->numVars(); i++) q[i] = q_trial[i];
554 if (missedOutput.empty())
return sigma;
560 if (!missedOutput.empty())
563 for (; iter != missedOutput.
end(); iter++)
566 sys->output_func(q_trial,event,yb);
571 sys->postStep(q_trial);
573 sys->output_func(q_trial,event,yb);
581 delete [] q;
delete [] q_trial;
delete [] event;
582 delete event_finder;
delete solver;
delete sys;
596 void tentative_step()
601 double step_size = solver->
integrate(q_trial,time_event);
603 bool state_event_exists =
604 event_finder->
find_events(event,q,q_trial,solver,step_size);
606 sigma = std::min<double>(step_size,time_event);
607 event[sys->
numEvents()] = time_event <= sigma;
608 event_exists =
event[sys->
numEvents()] || state_event_exists;
void gc_output(Bag< X > &gb)
Do not override. Invokes the ode_system gc_output method as needed.
Definition: adevs_hybrid.h:577
int numEvents() const
Get the number of state events.
Definition: adevs_hybrid.h:54
virtual ~Hybrid()
Destructor deletes everything.
Definition: adevs_hybrid.h:579
void output_func(Bag< X > &yb)
Do not override. Invokes the ode_system output function as needed.
Definition: adevs_hybrid.h:558
Hybrid(ode_system< X > *sys, ode_solver< X > *solver, event_locator< X > *event_finder)
Definition: adevs_hybrid.h:448
void der_func(const double *q, double *dq)
Do not override.
Definition: adevs_hybrid.h:210
Definition: adevs_hybrid.h:440
virtual ~ode_solver()
Destructor.
Definition: adevs_hybrid.h:396
void delta_conf(const Bag< X > &xb)
Definition: adevs_hybrid.h:534
virtual void gc_output(Bag< X > &gb)=0
Garbage collection function. This works just like the Atomic gc_output method.
double getWorseError() const
Definition: adevs_hybrid.h:203
virtual ~dae_se1_system()
Destructor.
Definition: adevs_hybrid.h:185
virtual ~ode_system()
Destructor.
Definition: adevs_hybrid.h:90
Definition: adevs_hybrid.h:406
ode_system< X > * getSystem()
Get the system that this solver is operating on.
Definition: adevs_hybrid.h:466
virtual void der_func(const double *q, double *dq)=0
Compute the derivative for state q and put it in dq.
void init(double *q)
Do not override.
Definition: adevs_hybrid.h:205
virtual double integrate(double *q, double h_lim)=0
virtual void postTrialStep(double *q, double *a)
Definition: adevs_hybrid.h:171
iterator begin() const
Get an iterator pointing to the first element in the bag.
Definition: adevs_bag.h:128
virtual void external_event(double *q, double e, const Bag< X > &xb)=0
The external transition function.
void confluent_event(double *q, const bool *state_event, const Bag< X > &xb)
Do not override.
Definition: adevs_hybrid.h:252
void internal_event(double *q, const bool *state_event)
Do not override.
Definition: adevs_hybrid.h:234
virtual void postTrialStep(double *q)
Definition: adevs_hybrid.h:74
void solve(const double *q)
Definition: adevs_hybrid.h:291
bool eventHappened() const
Did a discrete event occur at the last state transition?
Definition: adevs_hybrid.h:468
iterator end() const
Get an iterator to the end of the bag (i.e., just after the last element)
Definition: adevs_bag.h:130
dae_se1_system(int N_vars, int M_event_funcs, int A_alg_vars, double err_tol=1E-10, int max_iters=30, double alpha=-1.0)
Definition: adevs_hybrid.h:119
double getAlgVar(int i) const
Get the algebraic variables.
Definition: adevs_hybrid.h:138
int numVars() const
Get the number of state variables.
Definition: adevs_hybrid.h:52
void delta_ext(T e, const Bag< X > &xb)
Definition: adevs_hybrid.h:496
ode_solver(ode_system< X > *sys)
Definition: adevs_hybrid.h:384
Definition: adevs_fmi.h:57
virtual void postStep(double *q)
Definition: adevs_hybrid.h:68
double time_event_func(const double *q)
Override only if you have no time events.
Definition: adevs_hybrid.h:222
const double * getState() const
Get the array of state variables.
Definition: adevs_hybrid.h:464
virtual void internal_event(double *q, const bool *state_event)=0
The internal transition function.
void state_event_func(const double *q, double *z)
Override only if you have no state event functions.
Definition: adevs_hybrid.h:216
void delta_int()
Definition: adevs_hybrid.h:473
void postStep(double *q)
Do not override.
Definition: adevs_hybrid.h:228
virtual ~event_locator()
Destructor.
Definition: adevs_hybrid.h:427
ode_system(int N_vars, int M_event_funcs)
Make a system with N state variables and M state event functions.
Definition: adevs_hybrid.h:49
void external_event(double *q, double e, const Bag< X > &xb)
Do not override.
Definition: adevs_hybrid.h:243
virtual void state_event_func(const double *q, double *z)=0
Compute the state event functions for state q and put them in z.
Definition: adevs_hybrid.h:108
double getState(int k) const
Get the value of the kth continuous state variable.
Definition: adevs_hybrid.h:462
virtual double time_event_func(const double *q)=0
Compute the time event function using state q.
int getIterFailCount() const
Definition: adevs_hybrid.h:197
event_locator(ode_system< X > *sys)
Definition: adevs_hybrid.h:413
void insert(const T &t)
Put t into the bag.
Definition: adevs_bag.h:153
int numAlgVars() const
Get the number of algebraic variables.
Definition: adevs_hybrid.h:136
Definition: adevs_models.h:47
T ta()
Do not override.
Definition: adevs_hybrid.h:552
virtual bool find_events(bool *events, const double *qstart, double *qend, ode_solver< X > *solver, double &h)=0
virtual void init(double *q)=0
Copy the initial state of the model to q.
Definition: adevs_hybrid.h:45
virtual void confluent_event(double *q, const bool *state_event, const Bag< X > &xb)=0
The confluent transition function.
Definition: adevs_hybrid.h:377
void output_func(const double *q, const bool *state_event, Bag< X > &yb)
Do not override.
Definition: adevs_hybrid.h:261
virtual void output_func(const double *q, const bool *state_event, Bag< X > &yb)=0
The output function.