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