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: 2253
Valid CSS Valid XHTML 1.0 Strict