/[ascend]/trunk/models/steam/stab.py
ViewVC logotype

Contents of /trunk/models/steam/stab.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1328 - (show annotations) (download) (as text)
Thu Mar 8 03:05:27 2007 UTC (12 years, 11 months ago) by jpye
File MIME type: text/x-python
File size: 2637 byte(s)
Added stability analysis example in models/steam.
'unattached' icon was created some time ago but is not yet in use
1 # this is a script that computes the stability of the DAE system of equations
2 # by using the sparse matrix routines in scipy and plotting with matplotlib.
3 #
4 # you could get fancy and produce a root locus using this technique...
5
6 import ascpy
7
8 L = ascpy.Library()
9 L.load('steam/dsgsat3.a4c')
10 T = L.findType('dsgsat3')
11 M = T.getSimulation('sim',False)
12 M.run(T.getMethod('on_load'))
13 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
14 M.run(T.getMethod('configure_dynamic'))
15 M.solve(ascpy.Solver('QRSlv'),ascpy.SolverReporter())
16 T = L.findType('dsgsat3')
17 M.run(T.getMethod('free_states'))
18 # here is the peturbation...
19 M.qdot_s.setRealValueWithUnits(6000,"W/m")
20 # IDA has its own initial conditions solver, so no need to call QRSlv here
21 I = ascpy.Integrator(M)
22 I.setEngine('IDA')
23 I.setParameter('linsolver','DENSE')
24 I.setParameter('safeeval',True)
25 I.setParameter('rtol',1e-4)
26 I.setParameter('atolvect',False)
27 I.setParameter('atol',1e-4)
28 I.setParameter('maxord',3)
29 I.setInitialSubStep(0.001)
30 I.setReporter(ascpy.IntegratorReporterConsole(I))
31 I.setLogTimesteps(ascpy.Units("s"), 0.001, 3600, 10)
32 I.analyse()
33 F = file('ga.mm','w')
34 I.writeMatrix(F,'dg/dya')
35 F = file('gd.mm','w')
36 I.writeMatrix(F,'dg/dyd')
37 F = file('fa.mm','w')
38 I.writeMatrix(F,'df/dya')
39 F = file('fd.mm','w')
40 I.writeMatrix(F,'df/dyd')
41 F = file('fdp.mm','w')
42 I.writeMatrix(F,'df/dydp')
43 #I.solve()
44
45 from scipy import io
46 from scipy import linalg
47
48 gz = io.mmread('ga.mm')
49 gx = io.mmread('gd.mm')
50 fz = io.mmread('fa.mm')
51 fx = io.mmread('fd.mm')
52 fxp = io.mmread('fdp.mm')
53
54 print "gz", gz.shape
55 print "gx", gx.shape
56 print "fz", fz.shape
57 print "fx", fx.shape
58 print "fxp", fxp.shape
59
60 #import pylab
61
62 # dg/dy_a
63
64 #pylab.spy2(ga.todense())
65 #pylab.title("${dg}/{dy_a}$")
66 #pylab.show()
67
68 invgz = linalg.inv(gz.todense())
69
70 #pylab.figure()
71 #pylab.spy(invgz)
72 #pylab.title("$({dg}/{dy_d})^{-1}$")
73 #pylab.show()
74
75 # dg/dy_d
76
77 #pylab.figure()
78 #pylab.spy2(gd.todense())
79 #pylab.title("${dg}/{dy_d}$")
80 #pylab.show()
81
82 # df/dyd'
83
84 #pylab.figure()
85 #pylab.spy2(fdp.todense())
86 #pylab.title("${df}/{d\dot{y}_d}$")
87 #pylab.show()
88
89 invfxp = linalg.inv(fxp.todense())
90
91 #pylab.spy2(invfdp)
92 #pylab.title("$({df}/{dy_dp})^{-1}$")
93 #pylab.show()
94
95 dya_dyd = invgz * gx
96
97 print "gz^-1 gx",dya_dyd.shape
98
99 #pylab.spy2(dya_dyd.todense())
100 #pylab.title("${dy_a}/{dy_d}$")
101 #pylab.show()
102
103 B = fz * invgz * gx
104
105 print "fz gz^1 gz",B.shape
106 #pylab.spy2(fad.todense())
107 #pylab.title("${df}/{dy_a} * {dy_a}/{dy_d}$")
108 #pylab.show()
109
110 C = fx + B
111
112 D = - invfxp * C
113
114 e,v = linalg.eig(D.todense())
115
116 #print e
117
118 print "max re(e)",max(e.real)
119 print "min re(e)",min(e.real)
120
121 print "max im(e)",max(e.imag)
122 print "min in(e)",min(e.imag)
123
124 I.solve()
125

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