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;
71 const bool* state_event) = 0;
79 virtual void output_func(
const double *q,
const bool* state_event,
114 double err_tol = 1E-10,
int max_iters = 30,
double alpha = -1.0):
117 max_iters(max_iters),
124 atmp =
new double[A];
126 f[1] =
new double[A];
127 f[0] =
new double[A];
137 virtual void init(
double* q,
double* a) = 0;
144 virtual void alg_func(
const double* q,
const double* a,
double* af) = 0;
149 virtual void der_func(
const double* q,
const double* a,
double* dq) = 0;
153 virtual void state_event_func(
const double* q,
const double* a,
double* z) = 0;
161 virtual void postStep(
double* q,
double* a) = 0;
164 const bool* state_event) = 0;
167 double e,
const Bag<X>& xb) = 0;
170 const bool* state_event,
const Bag<X>& xb) = 0;
172 virtual void output_func(
const double *q,
const double* a,
173 const bool* state_event,
Bag<X>& yb) = 0;
264 void solve(
const double* q);
266 const int A, max_iters;
267 const double err_tol, alpha;
280 template <
typename X>
283 int iter_count = 0, alt, good;
284 double prev_err, err = 0.0, ee, beta, g2, alpha_tmp = alpha;
290 _adevs_dae_se_1_system_solve_try_it_again:
295 alg_func(q,a,f[alt]);
296 for (
int i = 0; i < A; i++)
304 a[i] += alpha_tmp*d[i];
308 while (iter_count < max_iters)
313 alg_func(q,a,f[good]);
315 for (
int i = 0; i < A; i++)
320 ee = fabs(f[good][i]);
321 if (ee > err) err = ee;
324 if (err < err_tol)
return;
329 for (
int i = 0; i < A; i++)
332 if (alpha_tmp < 0.0) alpha_tmp = -alpha_tmp;
333 else alpha_tmp *= -0.5;
334 goto _adevs_dae_se_1_system_solve_try_it_again;
340 for (
int i = 0; i < A; i++)
341 g2 += f[alt][i]*f[alt][i];
342 for (
int i = 0; i < A; i++)
343 beta += f[good][i]*(f[good][i]-f[alt][i]);
346 for (
int i = 0; i < A; i++)
348 d[i] = beta*d[i]-f[good][i];
350 a[i] += alpha_tmp*d[i];
380 virtual double integrate(
double* q,
double h_lim) = 0;
384 virtual void advance(
double* q,
double h) = 0;
414 virtual bool find_events(
bool* events,
const double *qstart,
430 template <
typename X,
class T =
double>
class Hybrid:
440 sys(sys),solver(solver),event_finder(event_finder),
443 q =
new double[sys->
numVars()];
444 q_trial =
new double[sys->
numVars()];
446 event_exists =
false;
448 for (
int i = 0; i < sys->
numVars(); i++) q[i] = q_trial[i];
465 if (!missedOutput.
empty())
467 missedOutput.
clear();
472 event_happened = event_exists;
475 sys->internal_event(q_trial,event);
479 for (
int i = 0; i < sys->numVars(); i++) q[i] = q_trial[i];
488 bool state_event_exists =
false;
489 event_happened =
true;
493 for (
int i = 0; i < sys->numVars(); i++)
495 solver->advance(q_trial,e);
497 event_finder->find_events(event,q,q_trial,solver,e);
499 if (state_event_exists)
502 sys->confluent_event(q_trial,event,xb);
503 for (
int i = 0; i < sys->numVars(); i++)
507 if (!state_event_exists)
509 solver->advance(q,e);
513 sys->external_event(q,e+e_accum,xb);
517 for (
int i = 0; i < sys->numVars(); i++) q_trial[i] = q[i];
526 if (!missedOutput.
empty())
528 missedOutput.
clear();
529 if (sigma > 0.0) event_exists =
false;
532 event_happened =
true;
534 sys->confluent_event(q_trial,event,xb);
535 else sys->external_event(q_trial,e_accum+
ta(),xb);
538 for (
int i = 0; i < sys->numVars(); i++) q[i] = q_trial[i];
544 if (missedOutput.
empty())
return sigma;
550 if (!missedOutput.
empty())
553 for (; iter != missedOutput.
end(); iter++)
556 sys->output_func(q_trial,event,yb);
561 sys->postStep(q_trial);
563 sys->output_func(q_trial,event,yb);
571 delete [] q;
delete [] q_trial;
delete [] event;
572 delete event_finder;
delete solver;
delete sys;
586 void tentative_step()
589 double time_event = sys->time_event_func(q);
591 double step_size = solver->integrate(q_trial,time_event);
593 bool state_event_exists =
594 event_finder->find_events(event,q,q_trial,solver,step_size);
596 sigma = std::min(step_size,time_event);
597 event[sys->numEvents()] = time_event <= sigma;
598 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:567
virtual void der_func(const double *q, const double *a, double *dq)=0
int numEvents() const
Get the number of state events.
Definition: adevs_hybrid.h:54
bool empty() const
Same as size()==0.
Definition: adevs_bag.h:126
virtual ~Hybrid()
Destructor deletes everything.
Definition: adevs_hybrid.h:569
void output_func(Bag< X > &yb)
Do not override. Invokes the ode_system output function as needed.
Definition: adevs_hybrid.h:548
virtual void internal_event(double *q, double *a, const bool *state_event)=0
The internal transition function.
Hybrid(ode_system< X > *sys, ode_solver< X > *solver, event_locator< X > *event_finder)
Definition: adevs_hybrid.h:438
void der_func(const double *q, double *dq)
Do not override.
Definition: adevs_hybrid.h:200
Definition: adevs_hybrid.h:430
virtual void init(double *q, double *a)=0
virtual ~ode_solver()
Destructor.
Definition: adevs_hybrid.h:386
void delta_conf(const Bag< X > &xb)
Definition: adevs_hybrid.h:524
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:193
virtual ~dae_se1_system()
Destructor.
Definition: adevs_hybrid.h:175
virtual ~ode_system()
Destructor.
Definition: adevs_hybrid.h:84
Definition: adevs_hybrid.h:396
void clear()
Remove all of the elements from the bag.
Definition: adevs_bag.h:144
ode_system< X > * getSystem()
Get the system that this solver is operating on.
Definition: adevs_hybrid.h:456
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:195
virtual double integrate(double *q, double h_lim)=0
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:242
void internal_event(double *q, const bool *state_event)
Do not override.
Definition: adevs_hybrid.h:224
void solve(const double *q)
Definition: adevs_hybrid.h:281
bool eventHappened() const
Did a discrete event occur at the last state transition?
Definition: adevs_hybrid.h:458
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:113
double getAlgVar(int i) const
Get the algebraic variables.
Definition: adevs_hybrid.h:132
int numVars() const
Get the number of state variables.
Definition: adevs_hybrid.h:52
virtual double time_event_func(const double *q, const double *a)=0
Compute the time event function using state q and algebraic variables a.
void delta_ext(T e, const Bag< X > &xb)
Definition: adevs_hybrid.h:486
virtual void alg_func(const double *q, const double *a, double *af)=0
ode_solver(ode_system< X > *sys)
Definition: adevs_hybrid.h:374
virtual void external_event(double *q, double *a, double e, const Bag< X > &xb)=0
The external transition function.
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:212
const double * getState() const
Get the array of state variables.
Definition: adevs_hybrid.h:454
virtual void output_func(const double *q, const double *a, const bool *state_event, Bag< X > &yb)=0
The output function.
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:206
void delta_int()
Definition: adevs_hybrid.h:463
void postStep(double *q)
Do not override.
Definition: adevs_hybrid.h:218
virtual void advance(double *q, double h)=0
virtual ~event_locator()
Destructor.
Definition: adevs_hybrid.h:417
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
virtual void state_event_func(const double *q, const double *a, double *z)=0
virtual void confluent_event(double *q, double *a, const bool *state_event, const Bag< X > &xb)=0
The confluent transition function.
void external_event(double *q, double e, const Bag< X > &xb)
Do not override.
Definition: adevs_hybrid.h:233
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:102
double getState(int k) const
Get the value of the kth continuous state variable.
Definition: adevs_hybrid.h:452
virtual double time_event_func(const double *q)=0
Compute the time event function using state q.
int getIterFailCount() const
Definition: adevs_hybrid.h:187
event_locator(ode_system< X > *sys)
Definition: adevs_hybrid.h:403
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:130
Definition: adevs_models.h:47
virtual void postStep(double *q, double *a)=0
T ta()
Do not override.
Definition: adevs_hybrid.h:542
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:367
void output_func(const double *q, const bool *state_event, Bag< X > &yb)
Do not override.
Definition: adevs_hybrid.h:251
virtual void output_func(const double *q, const bool *state_event, Bag< X > &yb)=0
The output function.