/[ascend]/trunk/models/westerberg/ivpDAE/ivpStepN.a4c
ViewVC logotype

Contents of /trunk/models/westerberg/ivpDAE/ivpStepN.a4c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 889 - (show annotations) (download) (as text)
Wed Oct 11 21:01:41 2006 UTC (13 years, 9 months ago) by aw0a
File MIME type: text/x-ascend
File size: 19434 byte(s)
adding simplified ivpDAE models directory and file
1 REQUIRE "atoms.a4l";
2 REQUIRE "plot.a4l";
3
4 (*
5 * ivpStep.a4c
6 * by Arthur W. Westerberg
7 * Part of the ASCEND Library
8 * $Date: 02/07/28 $
9 * $Revision: 1.7 $
10 * $Author: mthomas $
11 * $Source: /afs/cs.cmu.edu/project/ascend/Repository/models/ivpsystem.a4l,v $
12 *
13 * This file is part of the ASCEND Modeling Library.
14 *
15 * Copyright (C) 1994 - 1998 Carnegie Mellon University
16 *
17 * The ASCEND Modeling Library is free software; you can redistribute
18 * it and/or modify it under the terms of the GNU General Public
19 * License as published by the Free Software Foundation; either
20 * version 2 of the License, or (at your option) any later version.
21 *
22 * The ASCEND Modeling Library is distributed in hope that it
23 * will be useful, but WITHOUT ANY WARRANTY; without even the implied
24 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
25 * See the GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with the program; if not, write to the Free Software
29 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check
30 * the file named COPYING.
31 *)
32
33 (*============================================================================
34
35 I V P S T E P . A 4 C
36 -----------------------------
37
38 AUTHOR: Arthur W. Westerberg
39
40 DATES: 07/2002 - Original code (AWW)
41 07/2003 - Tested code. Works for simple test problem.
42 (AWW)
43 12/2003 - Added Tcl script that runs complete simulation
44 (AWW)
45
46 CONTENTS: Models for the numerical integration equations for
47 the two multistep methods: a bdf for stiff problems
48 and the Adams Moulton for non-stiff problems. There
49 is also a framework within which one creates the
50 model for the physical system. These models are for
51 taking one step of the independent variable.
52
53 ============================================================================*)
54
55
56 (* ---------------------------------------------------------- *)
57
58 MODEL ivpBase() REFINES cmumodel();
59
60 (* ----- ivpBase -----
61 The base model for solving DAE initial value problems using the
62 multistep methods: (a) a bdf method for stiff problems and (b)
63 the Adams Moulton method for non-stiff problems. A DAE ivp
64 model has two types of sets of equations:
65
66 (1) the model of the physical system at the k-th step
67
68 f(dydx(k), y(k), z(k), u(k), x(k)) = 0
69 h(y(k), z(k), u(k), x(k)) = 0
70 u(k) = u(x(k))
71
72 where y are the state, z the algebraic and u the forcing
73 variable vectors and x the scalar independent variable, and
74 where f is a vector of functions equal to the number of
75 states, h is a vector of algebraic equations equal in number
76 to the algebraic variables, and where one typically
77 specifies the forcing variables as functions only of the
78 independent variable.
79
80 (2) the numerical multistep integration equations
81
82 y(k) = y(k-1) + s(x[k]-x[k-1];
83 y[k-1],y[k-2],...;
84 dydx[k],dydx[k-1],...)
85
86 These equations relate the states at the current point to
87 the time derivatives at current and past points and states
88 at past points. They are not related to the model of the
89 physical system but only to the numerical method being used.
90
91 A common solving approach would be to solve the integration
92 equations in an outside loop, and, for each iteration of that
93 outside loop, use a Newton-based method to solve the physical
94 system model for the state derivatives and the algebraic
95 variables. Often this approach solves the outer loop using a
96 simple substitution method - which can be very slow to
97 converge. The following figure illustrates.
98
99
100 ------------------
101 u(x) ->---| (inner loop) |---- z -->
102 | solve model eqns |
103 ------>---| |-- dydx-->-
104 | ------------------ |
105 | |
106 | ------------------ |
107 | | (outer loop) | |
108 --- y,x =-| iterate inte- |---<-------
109 | gration eqns |
110 | |
111 ------------------
112
113
114 By creating both the model and integration equations through
115 ASCEND models, ASCEND can simultaneously converge the model
116 and integration equations using its sparse Newton-based solver
117 - thus making the solving much more efficient. We wrote this
118 model precisely to gain this efficiency.
119
120
121 The integration module - degrees of freedom considerations
122
123 For the q-th order bdf and Adams-Moulton (am) methods used
124 here, we write a q-th order polynomial in the independent
125 variable x for each state variable y[k][i]. This polynomial is
126 written in the form of a Taylor Series expanded around the
127 current point x[k].
128
129 A q-th order polynomial has the q+1 polynomial coefficients in
130 it a[k][i]/i!. The a[k][i] are estimate the q-th derivative of
131 the function at the current point, with the q-th term being an
132 estimate of the error in the fit. We determine these
133 coefficients by fitting this polynomial and/or the first
134 derivative of this polynomial to past y[j][i] and dy/dt[j][i],
135 j<k, data, respectively, and to the computed value for
136 dy/dt[k][i] at the current point. We would need to write
137 precisely q+1 such equations for this fitting. The bdf and am
138 methods differ in which equations we write to do this fitting.
139 See comments in later in the code for the actual polynomials we
140 write.
141
142 We then use this polynomial evaluated at the current point to
143 estimate the state y[k][i] at the current point x[k].
144
145 Looking back at the two boxes in the previous figure, the lower
146 integration box will have these q+1 equations in it, tabulated
147 values for previous points, and the means to use past data and
148 the current estimate for dy/dt[k][i] (the input we show above
149 for this box) to fit the coefficients for each state variable.
150 It will also write the polynomial one additional time to
151 provide the estimate for y[k][i] at the current point, which is
152 the output from this box. It thus contains q+2 polynomial
153 equations in it, in terms of q+1 unknown coefficients for each
154 state. The model equations introduce the states y[k][i] and their
155 derivatives dy/dx[k][i] for the current point and the equation to
156 compute the latter in terms of the former, often in the form
157
158 dy/dt[k][i] = f(all y[k],t).
159
160 So it introduces, in principle two variables and one equation
161 for each state.
162
163 To summarize, the integration box introduces q+1 equations and
164 q+2 variables per state while the model box introduces 1
165 equation and 2 new variables. The model is in square and
166 solvable.
167
168
169 Predictor behavior
170
171 The polynomial we fit for the current point can easily be used
172 to predict y[k+1][i], the values for the states i at the next
173 time step. We can also create polynomials to predict any of
174 the non state variables in the problem. We can include these
175 non state variables in the integration model and use the
176 current and q past values to determine the coefficients for a
177 Taylor Series polynomial that passes through these points.
178 This polynomial would only be useful for predicting as it play
179 no role in solving the problem. But often predicting one of
180 the algebraic variables might improve the chance to converge
181 the equations.
182
183
184 Error control
185
186 The package checks periodically to see discover the size step
187 it should take to keep the error at that the user has
188 prescribed. It can also change the order of the polynomial up
189 or down by one if such a move will allow for a larger step.
190
191 We write the polynomial in the form of a Taylor Series around
192 the current point. So the error in estimating the current
193 point, as a function of the step size we take, is the highest
194 order term in that series. We back compute the step size
195 needed to obtain the error we allow and use it as the step size
196 for the next integration step. We keep this step size for the
197 next q steps before checking again. Changing the order of the
198 polynomial can also affect the error. This package will change
199 the order up or down by one (i.e., to q-1 or q+1 (keeping in
200 the limits of 1 and 6) if that change allows it to take a
201 larger step. The next to last term in the Taylor Series
202 estimates the error when using a polynomial of order q-1. We
203 estimate the derivative for order q+1 as the difference in the
204 q-th derivative between the current point and the previous
205 point divided by the change in x between these two points -
206 i.e., a divided difference which is essentially a free
207 calculation.
208
209
210 Stopping conditions
211
212 We have created this package so one can stop whenever any of
213 the variables we are predicting (all the states and whichever
214 of the non states we select) crosses through zero. In taking a
215 normal step in this package, we use the error control criterion
216 to fix the step size in the independent variable, x, prior to
217 taking the next step. If after we solve the model for this
218 step, any one of the variables we have selected to use as
219 stopping conditions changes sign, we fix its value to zero and
220 resolve the entire model to compute the step size. As more
221 than one such variable may change sign in this final interval,
222 we test after each solution to see if another variable changed
223 sign in the reduced interval we just computed. If so, we
224 repeat by fixing this new variable to zero and resolving. In
225 this manner we will find the variable that changed sign at the
226 earliest point in the final interval.
227
228 *)
229
230 (* ----- ivpBase -----
231 TO USE THIS PACKAGE, the user must provide three models that
232 are refinements of models we provide here. The first two
233 will be refinements of integrationPoint, while the third must be
234 a refinement of ivpStep. See the comments for those models.
235 Also see the model ivpStep.testModel.a4c, which is a prototype
236 for what the user must write.
237
238 NOTE: THE MODELS IN THIS LIBRARY ARE WRITTEN USING DIMENSIONLESS
239 QUANTITIES. THE USER MUST CONVERT ANY DIMENSIONAL STATES
240 AND PREDICTED ALGEBRAIC VARIABLES, AND ALL STOPPING CONDITIONS,
241 ERRORS, ETC., TO DIMENSIONLESS VARIABLES FOR TYING INTO THESE
242 MODELS.
243
244 *)
245
246
247 END ivpBase;
248
249 (* ---------------------------------------------------------- *)
250
251 MODEL factorial(
252 N IS_A integer_constant;
253 )
254 WHERE(
255 N>=0;
256 N<=12; (* Note that 12! = 0.48E9, 2^31=2.1E9 *)
257 ) REFINES ivpBase();
258
259 (* create a vector of factorials, where fac[i] contains
260 i!, where i goes from 0 to N *)
261
262 fac[0..N] IS_A factor;
263
264 METHODS
265
266 (* ----- factorial ----- *)
267
268 METHOD default_self;
269 fac[0] := 1;
270
271 (* if N=0, the following loop produces no code *)
272 FOR i IN [1..N] DO
273 fac[i] := i*fac[i-1];
274 END FOR;
275 END default_self;
276
277 (* ----- factorial ----- *)
278
279 METHOD specify;
280 FIX fac[0..N];
281 END specify;
282
283 END factorial;
284
285 (* ---------------------------------------------------------- *)
286
287 MODEL pTrow(
288 set IS_A set OF integer_constant;
289 ) REFINES ivpBase();
290
291 (* This model is needed and instanced in
292 "MODEL pascalTriangle." No user should
293 instance it by itself. *)
294
295 col[set] IS_A factor;
296
297 END pTrow;
298
299 (* ---------------------------------------------------------- *)
300
301 MODEL pascalTriangle(
302 dim IS_A integer_constant;
303 ) REFINES ivpBase();
304
305 (* Creates a Pascal Triangle of dimension "dim."
306 A Pascal triangle gives the coefficients for the
307 polynomials resulting from expanding (x-y)^q,
308 for q = 0..dim. For dim=5, the triangle we form
309 here is
310
311 row(1): 1 1 1 1 1 1
312 row(2): 1 2 3 4 5
313 row(3): 1 3 6 10
314 row(4): 1 4 10
315 row(5): 1 5
316 row(6): 1
317
318 The coefficients for q=4 for (x-y)^4 are then in
319 the next to last place in each row, namely
320 1, 4, 6, 4, 1
321 as in
322 (1)x^4 - (4)x^3y + (6)x^2y^2 - (4)xy^3 + (1)y^4.
323 *)
324
325 FOR i IN [0..dim] CREATE
326 row[i] IS_A pTrow([0..dim-i]);
327 END FOR;
328
329 METHODS
330
331 (* ----- pascalTriangle ----- *)
332
333 METHOD default_self;
334 FOR i IN [0..dim] DO
335 row[0].col[i] := 1;
336 row[i].col[0] := 1;
337 END FOR;
338 FOR i IN [1..dim] DO
339 FOR j IN [1..dim-i] DO
340 row[i].col[j] := row[i-1].col[j]+row[i].col[j-1];
341 END FOR;
342 END FOR;
343 END default_self;
344
345 (* ----- pascalTriangle ----- *)
346
347 METHOD specify;
348 FOR i IN [0..dim] DO
349 FOR j IN [0..dim-i] DO
350 FIX row[i].col[j];
351 END FOR;
352 END FOR;
353 END specify;
354
355 END pascalTriangle;
356
357 (* ---------------------------------------------------------- *)
358 (* ---------------------------------------------------------- *)
359 (* ---------------------------------------------------------- *)
360 (* ---------------------------------------------------------- *)
361 (* ---------------------------------------------------------- *)
362 (* ---------------------------------------------------------- *)
363
364 MODEL yPolysDef(
365 y WILL_BE factor;
366 dx WILL_BE factor;
367 dim WILL_BE integer_constant;
368 a[0..dim] WILL_BE factor;
369 ) WHERE(
370 (dim >= 3) == TRUE;
371 );
372
373 (*
374 create a polynomial equation for y of order dim at the
375 independent variable dx.
376 *)
377
378 f IS_A factorial(dim);
379 yTerm[0..dim] IS_A factor;
380
381 yTerm[dim] = a[dim]/f.fac[dim];
382 FOR i IN [2..dim] CREATE
383 yTermEqn[i-1]: yTerm[i-1] =
384 a[i-1]/f.fac[i-1] + dx*yTerm[i];
385 END FOR;
386 yTermEqn0: y = a[0] + dx*yTerm[1];
387
388 END yPolysDef;
389
390 (* ---------------------------------------------------------- *)
391
392 MODEL dydxPolysDef(
393 dydx WILL_BE factor;
394 dx WILL_BE factor;
395 dim WILL_BE integer_constant;
396 a[0..dim] WILL_BE factor;
397 ) WHERE(
398 (dim >= 3) == TRUE;
399 );
400
401 (*
402 create a polynomial equation for dydx of order dim at the
403 independent variable dx.
404 *)
405
406 f IS_A factorial(dim);
407 dydxTerm[0..dim-1] IS_A factor;
408
409 dydxTerm[dim] = a[dim]/f.fac[dim-1];
410 FOR i IN [3..dim] CREATE
411 dydxTermEqn[i-1]: dydxTerm[i-1] =
412 a[i-1]/f.fac[i-2] + dx*dydxTerm[i];
413 END FOR;
414 dydxTermEqn1: dydx = a[1] + dx*dydxTerm[2];
415
416 END dydxPolysDef;
417
418 (* ---------------------------------------------------------- *)
419
420 MODEL state(
421 y WILL_BE solver_var;
422 dydx WILL_BE solver_var;
423 dim WILL_BE integer_constant;
424 ) WHERE(
425 (dim >= 3) == TRUE;
426 );
427
428 (*
429 Model to provide storage for the state variable y, its derivative
430 dydt, and the independent variable dx, written at current and
431 past points. Note that dx = 0 at the current point. All variables
432 are stored in dimensionless form (of type factor), by dividing
433 values by nominal values when storing or retrieving.
434
435 It also provides polynomials of order dim at current and past points
436 for the integration equations written "around" the current point.
437 *)
438
439
440 (* storage *)
441 yStore[0..dim] IS_A factor;
442 dydxStore[0..dim] IS_A factor;
443 dxStore[0..dim] IS_A factor;
444
445 (* polynomials *)
446 a[0..dim] IS_A factor;
447 FOR i IN [1..dim] CREATE
448 yPolys[i] IS_A yPolysDef(yStore[i], dxStore[i], dim, a);
449 dydxPolys[i-1] IS_A yPolysDef(dydxStore[i-1], dxStore[i-1], dim, a);
450 END FOR;
451
452 (* Note that a[0] and a[1] predict y and dydx at the current point *)
453 yPredEqn: y/y.nominal = a[0];
454 dydxPredEqn: dydx/dydx.nominal = a[1];
455
456 (*
457
458 The solver can now choose which order and which method to use as
459 follows for an order q method
460
461 Set coefficients a[q+1..dim] to zero and fix.
462
463 To choose the AM method, set only the following included flags to TRUE
464 polys.dydxTermEqn[0..q-1]
465 polys.yTermEqn[1]
466
467 To choose the BDF method, set only the following included flags to TRUE
468 polys.dydxTermEqn[0]
469 polys.yTermEqn[1..q]
470
471 One will have written q+1 eqns in either case in the q+1 coefficients
472 a[0..q].
473
474 *)
475
476 METHODS
477
478 METHOD default_self;
479 END default_self;
480
481 METHOD specify;
482 END specify;
483
484 METHOD values;
485 END values;
486
487 METHOD setUseMethodAM;
488 END setUseMethodAM;
489
490 METHOD setUseMethodBDF;
491 END setUseMethodBDF;
492
493 METHOD incrementPolyOrder;
494 END incrementPolyOrder;
495
496 METHOD decrementPolyOrder;
497 END decrementPolyOrder;
498
499 METHOD stepX;
500 END stepX;
501
502 END state;
503
504 (* ---------------------------------------------------------- *)
505
506 MODEL predictedAlgebraic(
507 z WILL_BE solver_var;
508 dim WILL_BE integer_constant;
509 ) WHERE(
510 (dim >= 3) == TRUE;
511 );
512
513 (* storage *)
514 zStore[0..dim] IS_A factor;
515 dxStore[0..dim] IS_A factor;
516
517 (* polynomials *)
518 a[0..dim] IS_A factor;
519 FOR i IN [0..dim] CREATE
520 zPolys[i] IS_A yPolysDef(zStore[i], dxStore[i], dim, a);
521 END FOR;
522
523 (*
524 Note that algebraic prediction is only to improve initial guesses
525 for z at the current point so z will not necessarily equal its
526 converged current value. Do not include a prediction equation.
527
528 The solver can now choose which order to use as follows for an order
529 q method
530
531 Set coefficients a[q+1..dim] to zero and fix.
532
533 Set only the following include flags to TRUE
534 polys.yTermEqn[0..q]
535
536 One will have written q+1 eqns in the q+1 coefficients a[0..q].
537 *)
538
539 END predictedAlgebraic;
540
541 (* ---------------------------------------------------------- *)
542
543 MODEL myDynamicModel;
544
545 maxPolyOrder IS_A integer_constant;
546 maxPolyOrder :== 6;
547
548 t,dt IS_A time;
549 holdup IS_A mole;
550 dholdup_dt IS_A delta_molar_rate;
551 flowIn, flowOut,
552 flowStop, flowMin IS_A molar_rate;
553 Cv IS_A positive_factor; (* valve constant *)
554
555 holdupState IS_A state(holdup, dholdup_dt, maxPolyOrder);
556 predFlowStop IS_A predictedAlgebraic(flowStop, maxPolyOrder);
557
558 eqnDynTankholdup: dholdup_dt = flowIn - flowOut;
559 eqnFlowOut: (flowOut/1{mol/s})^2 = Cv^2*(holdup/1.0{mol});
560 eqnStoppingCond: flowStop = flowOut - flowMin;
561
562 METHODS
563
564 METHOD default_self;
565 END default_self;
566
567 METHOD specifyForInitializing;
568 END specifyForInitializing;
569
570 METHOD valuesForInitializing;
571 END valuesForInitializing;
572
573 METHOD specifyForStepping;
574 END specifyForStepping;
575
576 METHOD valuesForStepping;
577 END valuesForStepping;
578
579 METHOD testForIndexProblem;
580 END testForIndexProblem;
581
582 END myDynamicModel;
583
584 (* ---------------------------------------------------------- *)

Properties

Name Value
svn:executable *

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22