Parent Directory | 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)

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 | (* ---------------------------------------------------------- *) |

Name | Value |
---|---|

svn:executable |
* |

john.pye@anu.edu.au | ViewVC Help |

Powered by ViewVC 1.1.22 |