```   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)}
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
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}
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.}
100: writeln(output,'h (Seconds): ',h:12:7);
101:
102: {Get desired number of time steps for simulation.}
104: writeln(output,'maxsteps: ',maxsteps:4);
105:
106: {Get parameters for source V1 and initialize the source.}
107: with v1 do
108:    begin
110:    writeln(output,'V1 (Volts RMS): ',rtemp:9:3);
111:    ampl := rtemp * SQRT2;
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
123:    writeln(output,'R1 (Ohms): ',r:9:3);
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
139:    writeln(output,'R2 (Ohms): ',r:9:3);
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
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:
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;
354: for i := 1 to N do
355:    begin
357:    for j := 1 to N do
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
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.
```
 Last modified: 1984-04-01 Generated: 2016-12-26 Generated by src2html V0.67 page hit count: 544