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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1331 - (show annotations) (download) (as text)
Thu Mar 8 06:03:43 2007 UTC (12 years, 9 months ago) by jpye
File MIME type: text/x-python
File size: 2863 byte(s)
Reverted shm.a4c.
Added shmroots.a4c which displays the two complex roots of the SHM system using the extpy 'roots' module.
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.setParameter('calcic','YA_YDP')
30 I.setInitialSubStep(0.001)
31 I.setReporter(ascpy.IntegratorReporterConsole(I))
32 I.setLogTimesteps(ascpy.Units("s"), 0.001, 0.002, 10)
33 I.analyse()
34 F = file('ga.mm','w')
35 I.writeMatrix(F,'dg/dya')
36 F = file('gd.mm','w')
37 I.writeMatrix(F,'dg/dyd')
38 F = file('fa.mm','w')
39 I.writeMatrix(F,'df/dya')
40 F = file('fd.mm','w')
41 I.writeMatrix(F,'df/dyd')
42 F = file('fdp.mm','w')
43 I.writeMatrix(F,'df/dydp')
44 #I.solve()
45
46 from scipy import io
47 from scipy import linalg
48
49 gz = io.mmread('ga.mm')
50 gx = io.mmread('gd.mm')
51 fz = io.mmread('fa.mm')
52 fx = io.mmread('fd.mm')
53 fxp = io.mmread('fdp.mm')
54
55 print "gz", gz.shape
56 print "gx", gx.shape
57 print "fz", fz.shape
58 print "fx", fx.shape
59 print "fxp", fxp.shape
60
61 #import pylab
62
63 # dg/dy_a
64
65 #pylab.spy2(ga.todense())
66 #pylab.title("${dg}/{dy_a}$")
67 #pylab.show()
68
69 invgz = linalg.inv(gz.todense())
70
71 #pylab.figure()
72 #pylab.spy(invgz)
73 #pylab.title("$({dg}/{dy_d})^{-1}$")
74 #pylab.show()
75
76 # dg/dy_d
77
78 #pylab.figure()
79 #pylab.spy2(gd.todense())
80 #pylab.title("${dg}/{dy_d}$")
81 #pylab.show()
82
83 # df/dyd'
84
85 #pylab.figure()
86 #pylab.spy2(fdp.todense())
87 #pylab.title("${df}/{d\dot{y}_d}$")
88 #pylab.show()
89
90 invfxp = linalg.inv(fxp.todense())
91
92 #pylab.spy2(invfdp)
93 #pylab.title("$({df}/{dy_dp})^{-1}$")
94 #pylab.show()
95
96 dya_dyd = invgz * gx
97
98 print "gz^-1 gx",dya_dyd.shape
99
100 #pylab.spy2(dya_dyd.todense())
101 #pylab.title("${dy_a}/{dy_d}$")
102 #pylab.show()
103
104 B = fz * invgz * gx
105
106 print "fz gz^1 gz",B.shape
107 #pylab.spy2(fad.todense())
108 #pylab.title("${df}/{dy_a} * {dy_a}/{dy_d}$")
109 #pylab.show()
110
111 C = fx + B
112
113 D = - invfxp * C
114
115 e,v = linalg.eig(D.todense())
116
117 #print e
118
119 print "max re(e)",max(e.real)
120 print "min re(e)",min(e.real)
121
122 print "max im(e)",max(e.imag)
123 print "min in(e)",min(e.imag)
124
125 I.solve()
126
127 import pylab, sys
128 sys.stderr.write("about to plot...")
129 pylab.plot(e.real,e.imag,'rx')
130 pylab.show()
131 sys.stderr.write("DONE\n")
132
133 I.setLogTimesteps(ascpy.Units("s"), 0.002, 3600, 10)
134 I.solve()
135

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