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.
|
Last modified: 1984-04-01
Generated: 2016-12-26 |
Generated by src2html V0.67
page hit count: 547 |
|