![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to parall.p CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSD / ucb / pascal / pdx / test|
Return to parall.p CVS log | Up to [CSRG BSD Unix] / 43BSD / ucb / pascal / pdx / test |
1.1 root 1: program Parall(input,output);
2:
3: {Declare constants for unit conversions, convergence tests, etc.}
4:
5: const SQRT2 = 1.4142136; {Square root of 2}
6: TWOPI = 6.2831853; {Two times pi}
7: MINALPHA = 0.001; {Minimum y-step size}
8: ROUGHLYZERO = 0.001; {Approximation to zero for convergence}
9: YTHRESHOLD = 40.0; {Heuristic constant for thresholding}
10: N = 8; {Vector and matrix dimension}
11:
12:
13: {Declare types for circuit elements, vectors, and matrices.}
14:
15: type VSOURCE = record
16: ampl: real; {Volts (peak)}
17: freq: real; {Radians/second}
18: xindex: integer; {Index for x value}
19: yindex: integer; {Index for y value}
20: end;
21:
22: RLPAIR = record
23: r: real; {Ohms}
24: l: real; {Henries}
25: islope: real; {Amps/second}
26: invariant: real; {Trapezoidal invariant}
27: lasttime: real; {Most recent time}
28: xindex: array [1..2] of integer; {x value indices}
29: yindex: array [1..2] of integer; {y value indices}
30: end;
31:
32: CAPACITOR = record
33: c: real; {Farads}
34: vslope: real; {Volts/second}
35: invariant: real; {Trapezoidal invariant}
36: lasttime: real; {Most recent time}
37: xindex: array [1..2] of integer; {x value indices}
38: yindex: array [1..2] of integer; {y value indices}
39: end;
40:
41: VECTOR = array [1..N] of real;
42:
43: MATRIX = array [1..N,1..N] of real;
44:
45:
46: {Declare global variables for central simulation information.}
47:
48: var ground: VSOURCE; {Ground -- a fake source}
49: itcount: integer; {Main routine iteration count}
50: update: integer; {Update loop counter for main}
51: optimcount: integer; {Number of optimization steps}
52: newtoncount: integer; {Number of Newton steps}
53: v1: VSOURCE; {Voltage source V1}
54: rl1: RLPAIR; {R1/L1 resistor-inductor pair}
55: rl2: RLPAIR; {R2/L2 resistor-inductor pair}
56: c1: CAPACITOR; {Capacitor C1}
57: a: MATRIX; {Global matrix A}
58: b: MATRIX; {Global matrix B}
59: jac: MATRIX; {Global Jacobian matrix}
60: x: VECTOR; {Global vector of dependents}
61: xnext: VECTOR; {Next x-vector for simulation}
62: y: VECTOR; {Global vector of independents}
63: ynext: VECTOR; {Next y-vector for simulation}
64: gradient: VECTOR; {Global gradient vector}
65: h: real; {Time step value}
66: time: real; {Current time}
67: lastychange: real; {YStep's most recent y-change}
68: timestep: integer; {Current timestep number}
69: maxsteps: integer; {Number of time steps to run}
70: oldxnorm: real; {Old one-norm of x-vector}
71: newxnorm: real; {New one-norm of x-vector}
72: closenough: boolean; {Flag to indicate convergence}
73:
74:
75:
76:
77: {The following procedure initializes everything for the program based
78: on the little test circuit suggested by Sarosh. The user is asked
79: to specify the simulation and circuit parameters, and then the matrix
80: and vector values are set up.}
81:
82: procedure InitializeEverything;
83: var i,j: integer;
84: rtemp: real;
85: begin
86:
87: {Ready the input and output files (almost nil for Berkeley).}
88: writeln(output);
89: writeln(output,'*** Simulation Output Record ***');
90: writeln(output);
91: writeln(output);
92:
93: {Initialize convergence test/indication variables.}
94: oldxnorm := 0.0;
95: newxnorm := 0.0;
96: lastychange := 0.0;
97:
98: {Get desired time step size for simulation.}
99: readln(input,h);
100: writeln(output,'h (Seconds): ',h:12:7);
101:
102: {Get desired number of time steps for simulation.}
103: readln(input,maxsteps);
104: writeln(output,'maxsteps: ',maxsteps:4);
105:
106: {Get parameters for source V1 and initialize the source.}
107: with v1 do
108: begin
109: readln(input,rtemp);
110: writeln(output,'V1 (Volts RMS): ',rtemp:9:3);
111: ampl := rtemp * SQRT2;
112: readln(input,rtemp);
113: writeln(output,'f (Hertz): ',rtemp:9:3);
114: freq := rtemp * TWOPI;
115: xindex := 1;
116: yindex := 1;
117: end;
118:
119: {Get parameters for R1/L1 pair and initialize the pair.}
120: with rl1 do
121: begin
122: readln(input,r);
123: writeln(output,'R1 (Ohms): ',r:9:3);
124: readln(input,l);
125: writeln(output,'L1 (Henries): ',l:12:7);
126: islope := 0.0;
127: invariant := 0.0;
128: lasttime := -1.0; {Negative to force first update}
129: xindex[1] := 2;
130: yindex[1] := 2;
131: xindex[2] := 3;
132: yindex[2] := 3;
133: end;
134:
135: {Get parameters for R2/L2 pair and initialize the pair.}
136: with rl2 do
137: begin
138: readln(input,r);
139: writeln(output,'R2 (Ohms): ',r:9:3);
140: readln(input,l);
141: writeln(output,'L2 (Henries): ',l:12:7);
142: islope := 0.0;
143: invariant := 0.0;
144: lasttime := -1.0; {Negative to force first update}
145: xindex[1] := 4;
146: yindex[1] := 4;
147: xindex[2] := 5;
148: yindex[2] := 5;
149: end;
150:
151: {Get parameters for capacitor C1 and initialize the device.}
152: with c1 do
153: begin
154: readln(input,c);
155: writeln(output,'C1 (Farads): ',c:12:7);
156: vslope := 0.0;
157: invariant := 0.0;
158: lasttime := -1.0; {Negative to force first update}
159: xindex[1] := 6;
160: yindex[1] := 6;
161: xindex[2] := 7;
162: yindex[2] := 7;
163: end;
164:
165: {Initialize the ground node.}
166: with ground do
167: begin
168: ampl := 0.0;
169: freq := 0.0;
170: xindex := 8;
171: yindex := 8;
172: end;
173:
174: {Zero out all the vectors and matrices.}
175: for i := 1 to N do
176: begin
177: x[i] := 0.0;
178: y[i] := 0.0;
179: for j := 1 to N do
180: begin
181: a[i,j] := 0.0;
182: b[i,j] := 0.0;
183: jac[i,j] := 0.0;
184: end;
185: end;
186:
187: {Initialize the A matrix.}
188: a[1,2] := -1.0;
189: a[2,3] := 1.0;
190: a[2,4] := -1.0;
191: a[3,5] := 1.0;
192: a[4,7] := 1.0;
193: a[5,1] := 1.0;
194: a[7,6] := 1.0;
195: a[8,8] := 1.0;
196:
197: {Initialize the B matrix.}
198: b[1,1] := 1.0;
199: b[3,7] := -1.0;
200: b[4,8] := -1.0;
201: b[5,2] := 1.0;
202: b[6,3] := 1.0;
203: b[6,4] := 1.0;
204: b[7,5] := 1.0;
205: b[8,6] := 1.0;
206:
207: {Initialize the Jacobian matrix.}
208: rtemp := h / (2.0 * rl1.l + rl1.r * h);
209: jac[2,2] := rtemp;
210: jac[3,3] := rtemp;
211: jac[2,3] := -rtemp;
212: jac[3,2] := -rtemp;
213: rtemp := h / (2.0 * rl2.l + rl2.r * h);
214: jac[4,4] := rtemp;
215: jac[5,5] := rtemp;
216: jac[4,5] := -rtemp;
217: jac[5,4] := -rtemp;
218: jac[6,6] := -1.0;
219: jac[7,7] := 1.0;
220: jac[7,6] := h / (2.0 * c1.c);
221: end;
222:
223:
224:
225:
226: {The following procedure solves the equation Ax=b for an N x N system
227: of linear equations, where A is the coefficient matrix, b is the
228: right-hand-side vector, and x is the vector of unknowns. Gaussian
229: elimination with maximal column pivots is used. }
230:
231: procedure Solve(a: MATRIX; b: VECTOR; var x: VECTOR);
232: var y,z: real;
233: i,j,k,k1: integer;
234: begin
235: for k := 1 to N-1 do
236: begin
237: y := abs(a[k,k]);
238: j := k;
239: k1 := k + 1;
240: for i := k1 to N do
241: if abs(a[i,k]) > y then
242: begin
243: j := i;
244: y := abs(a[i,k]);
245: end;
246: for i := k to N do
247: begin
248: y := a[k,i];
249: a[k,i] := a[j,i];
250: a[j,i] := y;
251: end;
252: y := b[k];
253: b[k] := b[j];
254: b[j] := y;
255: z := a[k,k];
256: for i := k1 to N do
257: begin
258: y := a[i,k] / z;
259: a[i,k] := y;
260: for j := k1 to N do a[i,j] := a[i,j] - y * a[k,j];
261: b[i] := b[i] - y * b[k];
262: end;
263: end;
264: x[N] := b[N] / a[N,N];
265: for i := N-1 downto 1 do
266: begin
267: y := b[i];
268: for j := i+1 to N do y := y - a[i,j] * x[j];
269: x[i] := y / a[i,i];
270: end;
271: end;
272:
273:
274: {The following procedure computes and returns a vector called "deltay",
275: which is the change in the y-vector prescribed by the Newton-Rhapson
276: algorithm.}
277:
278: procedure NewtonStep(var deltay: VECTOR);
279: var phi: VECTOR;
280: m: MATRIX;
281: i,j,k: integer;
282: begin
283: for i := 1 to N do
284: begin
285: phi[i] := 0.0;
286: for j := 1 to N do
287: begin
288: phi[i] := phi[i] + a[i,j] * y[j] + b[i,j] * x[j];
289: m[i,j] := -a[i,j];
290: for k := 1 to N do m[i,j] := m[i,j] - b[i,k] * jac[k,j];
291: end;
292:
293: end;
294: Solve(m,phi,deltay);
295: end;
296:
297:
298:
299:
300: {The following function computes the value of theta, the objective
301: function, given the x and y vectors.}
302:
303: function ThetaValue(x,y: VECTOR): real;
304: var i,j: integer;
305: phielem: real;
306: theta: real;
307: begin
308: theta := 0.0;
309: for i:= 1 to N do
310: begin
311: phielem := 0.0;
312: for j := 1 to N do
313: phielem := phielem + a[i,j] * y[j] + b[i,j] * x[j];
314: theta := theta + phielem * phielem;
315: end;
316: ThetaValue := theta;
317: end;
318:
319:
320: {The following function computes the theta value associated with a
321: proposed step of size alpha in the direction of the gradient.}
322:
323: function Theta(alpha: real): real;
324: var ythere: VECTOR;
325: i: integer;
326: begin
327: for i := 1 to N do
328: ythere[i] := y[i] - alpha * gradient[i];
329: Theta := ThetaValue(x,ythere);
330: end;
331:
332:
333: {The following procedure computes the gradient of the objective
334: function (theta) with respect to the vector y.}
335:
336: procedure ComputeGradient;
337: var i,j,k: integer;
338: m: MATRIX;
339: v: VECTOR;
340: begin
341: {Compute v = Ay + Bx and M = A' + J'B'.}
342: for i := 1 to N do
343: begin
344: v[i] := 0.0;
345: for j := 1 to N do
346: begin
347: v[i] := v[i] + a[i,j] * y[j] + b[i,j] * x[j];
348: m[i,j] := a[j,i];
349: for k := 1 to N do
350: m[i,j] := m[i,j] + jac[k,i] * b[j,k];
351: end;
352: end;
353: {Compute gradient = 2Mv.}
354: for i := 1 to N do
355: begin
356: gradient[i] := 0.0;
357: for j := 1 to N do
358: gradient[i] := gradient[i] + m[i,j] * v[j];
359: gradient[i] := 2.0 * gradient[i];
360: end;
361: end;
362:
363:
364: {The following procedure computes the bounds on alpha, the step size
365: to take in the gradient direction. The bounding is done according
366: to an algorithm suggested in S.W.Director's text on circuits.}
367:
368: procedure GetAlphaBounds(var lower,upper: real);
369: var alpha: real;
370: oldtheta,newtheta: real;
371: begin
372: if ThetaValue(x,y) = 0.0
373: then
374: begin
375: lower := 0;
376:
377: upper := 0;
378: end
379: else
380: begin
381: lower := MINALPHA;
382: oldtheta := Theta(lower);
383: upper := MINALPHA * 2.0;
384: newtheta := Theta(upper);
385: if newtheta <= oldtheta then
386: begin
387: alpha := upper;
388: repeat
389: begin
390: oldtheta := newtheta;
391: alpha := alpha * 2.0;
392: newtheta := Theta(alpha);
393: end
394: until (newtheta > oldtheta);
395: lower := alpha / 4.0;
396: upper := alpha;
397: end;
398: end;
399: end;
400:
401:
402: {The following function computes the best value of alpha for the step
403: in the gradient direction. This best value is the value that minimizes
404: theta along the gradient-directed path.}
405:
406: function BestAlpha(lower,upper: real): real;
407: var alphaL,alphaU,alpha1,alpha2: real;
408: thetaL,thetaU,theta1,theta2: real;
409:
410: begin
411: alphaL := lower;
412: thetaL := Theta(alphaL);
413: alphaU := upper;
414: thetaU := Theta(alphaU);
415: alpha1 := 0.381966 * alphaU + 0.618034 * alphaL;
416: theta1 := Theta(alpha1);
417: alpha2 := 0.618034 * alphaU + 0.381966 * alphaL;
418: theta2 := Theta(alpha2);
419: repeat
420: if theta1 < theta2
421: then
422: begin
423: alphaU := alpha2;
424: thetaU := theta2;
425: alpha2 := alpha1;
426: theta2 := theta1;
427: alpha1 := 0.381966 * alphaU + 0.618034 * alphaL;
428: theta1 := Theta(alpha1);
429: end
430: else
431: begin
432: alphaL := alpha1;
433: thetaL := theta1;
434: alpha1 := alpha2;
435: theta1 := theta2;
436: alpha2 := 0.618034 * alphaU + 0.381966 * alphaL;
437: theta2 := Theta(alpha2);
438: end
439: until abs(thetaU - thetaL) <= ROUGHLYZERO;
440: BestAlpha := (alphaL + alphaU) / 2.0;
441: end;
442:
443:
444: {The following procedure computes and returns a vector called "deltay",
445: which is the change in the y-vector prescribed by the steepest-descent
446: approach.}
447:
448: procedure OptimizationStep(var deltay: VECTOR);
449: var lower,upper: real;
450: alpha: real;
451: i: integer;
452: begin
453: ComputeGradient;
454: GetAlphaBounds(lower,upper);
455: if lower <> upper then
456: begin
457: alpha := BestAlpha(lower,upper);
458: for i:= 1 to N do deltay[i] := - alpha * gradient[i];
459: end;
460: end;
461:
462:
463:
464:
465: {The following function computes the one-norm of a vector argument.
466: The length of the argument vector is assumed to be N.}
467:
468: function OneNorm(vec: VECTOR): real;
469: var sum: real;
470: i: integer;
471: begin
472: sum := 0;
473: for i := 1 to N do sum := sum + abs(vec[i]);
474: OneNorm := sum;
475: end;
476:
477:
478: {The following procedure takes a y-step, using the optimization
479: approach when far from the solution and the Newton-Rhapson approach
480: when fairly close to the solution.}
481:
482: procedure YStep;
483: var deltay: VECTOR;
484: ychange: real;
485: scale: real;
486: i: integer;
487: begin
488: NewtonStep(deltay);
489: ychange := OneNorm(deltay);
490: if ychange > YTHRESHOLD
491: then
492: {
493: begin
494: OptimizationStep(deltay);
495: ychange := OneNorm(deltay);
496: if ychange > YTHRESHOLD then
497: }
498: begin
499: scale := YTHRESHOLD/ychange;
500: for i := 1 to N do deltay[i] := scale * deltay[i];
501: optimcount := optimcount + 1;
502: end {;}
503: {
504: optimcount := optimcount + 1;
505: end
506: }
507: else
508: begin
509: newtoncount := newtoncount + 1;
510: end;
511: for i := 1 to N do ynext[i] := y[i] + deltay[i];
512: end;
513:
514:
515:
516:
517: {The following procedure updates the output of a voltage source
518: given the current time.}
519:
520: procedure VsourceStep(vn: VSOURCE);
521: begin
522: with vn do
523: xnext[xindex] := ampl * sin(freq * time);
524: end;
525:
526:
527: {The following procedure updates the outputs of a resistor-inductor
528: pair given the time step to take...that is, this routine takes a
529: time step for resistor-inductor pairs. The new outputs are found
530: by applying the trapezoidal rule.}
531:
532: procedure RLPairStep(var rln: RLPAIR);
533: begin
534: with rln do
535: begin
536: if (time > lasttime) then
537: begin
538: lasttime := time;
539: invariant := xnext[xindex[1]] + (h / 2.0) * islope;
540: end;
541: islope := (y[yindex[1]] - y[yindex[2]] - r * xnext[xindex[1]]) / l;
542: xnext[xindex[1]] := invariant + (h / 2.0) * islope;
543: xnext[xindex[2]] := - xnext[xindex[1]];
544: end;
545: end;
546:
547:
548: {The following procedure updates the outputs of a capacitor given the
549: time step...it takes the time step using a trapezoidal rule iteration.}
550:
551: procedure CapacitorStep(var cn: CAPACITOR);
552: var v: real;
553: begin
554: with cn do
555: begin
556: if (time > lasttime) then
557: begin
558: lasttime := time;
559: v := xnext[xindex[2]] - y[yindex[2]];
560: invariant := v + (h / 2.0) * vslope;
561: end;
562: vslope := y[yindex[1]] / c;
563: v := invariant + (h / 2.0) * vslope;
564: xnext[xindex[1]] := - y[yindex[1]];
565: xnext[xindex[2]] := y[yindex[2]] + v;
566: end;
567: end;
568:
569:
570: {The following procedure controls the overall x-step for the
571: specific circuit to be simulated.}
572:
573: procedure XStep;
574: begin
575: VsourceStep(v1);
576: RLPairStep(rl1);
577: RLPairStep(rl2);
578: CapacitorStep(c1);
579: VsourceStep(ground);
580: end;
581:
582:
583:
584:
585: {The following procedure prints out the values of all the voltages and
586: currents for the circuit at each time step.}
587:
588: procedure PrintReport;
589: begin
590: writeln(output);
591: writeln(output);
592: writeln(output,'REPORT at Time = ',time:8:5,' seconds');
593: writeln(output,'Number of iterations used: ',itcount:2);
594: writeln(output,'Number of truncations: ',optimcount:2);
595: writeln(output,'Number of Newton y-steps: ',newtoncount:2);
596: writeln(output,'The voltages and currents are:');
597: writeln(output,' V1 = ',x[1]:12:7,' I1 = ',y[1]:12:7);
598: writeln(output,' V2 = ',y[2]:12:7,' I2 = ',x[2]:12:7);
599: writeln(output,' V3 = ',y[3]:12:7,' I3 = ',x[3]:12:7);
600: writeln(output,' V4 = ',y[4]:12:7,' I4 = ',x[4]:12:7);
601: writeln(output,' V5 = ',y[5]:12:7,' I5 = ',x[5]:12:7);
602: writeln(output,' V6 = ',x[7]:12:7,' I6 = ',y[6]:12:7);
603: writeln(output,' V7 = ',y[7]:12:7,' I7 = ',x[6]:12:7);
604: writeln(output,' V8 = ',x[8]:12:7,' I8 = ',y[8]:12:7);
605: end;
606:
607:
608:
609:
610: {Finally, the main routine controls the whole simulation process.}
611:
612: begin
613: InitializeEverything;
614: PrintReport;
615: for timestep := 1 to maxsteps do
616: begin
617: itcount := 0;
618: optimcount := 0;
619: newtoncount := 0;
620: closenough := FALSE;
621: time := h * timestep;
622: repeat
623: begin
624: itcount := itcount + 1;
625: YStep;
626: XStep;
627: for update := 1 to N do
628: begin
629: x[update] := xnext[update];
630: y[update] := ynext[update];
631: end;
632: oldxnorm := newxnorm;
633: newxnorm := OneNorm(x);
634: if abs(newxnorm - oldxnorm) <= ROUGHLYZERO
635: then closenough := TRUE;
636: end;
637: if itcount < 4 then closenough := FALSE;
638: until (itcount = 25) or closenough;
639: PrintReport;
640: end;
641: end.
642:
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.