Appendix A

Meta Calculus Holonic Modeling

 

Western culture, with its hyper-individualism philosophy has been more vulnerable to the excesses of reductionism and diversity than Eastern culture whose religious roots are founded in holism and unity. The notion of balance, so fundamental in Eastern thought and so neglected in Western thought is lacking in Western science because reductionist math does not reach the level in Western education where balance is the overriding process of solution. When everything is reduced to a single equation in one unknown—the myopic funneling effect of our education pedagogy and the enforced discipline of our scientific modeling—the idea of balance becomes secondary. In effect it is assumed not to matter, because it does not appear in the equation, having been neglected in the reductionistic model-building process.

The computer has always held the promise to extricate us from this grip of reductionism by bringing us the means to solve holistic models—multidimensional, multivariable systems of equations wherein balance is the essence of solution. But its administration by the computer technocracy has delayed the onset of holistic modeling in education and in science.

In This Appendix

Holonic Theory

Meta Calculus Holonic Architecture

Basic Holonic Templates

Chemical Equilibrium – Nonlinear Equations

Admittance Circuit Fitting – Scientific Method

Nested Holonic Templates – Design Patterns

Pilot Ejection Altitude-Speed Profiling

Bipropellant Rocket – Implicit Differential Equations

Three-Stage Rocket Design Optimization

Chemical Kinetics – Scientific Method

Optimal Design and Control (ODC) Template

ODC Template Expanded – Wing Design Optimization

Systems Engineering End-User Perspective

Meta Calculus in Education

 

Holonic Theory

Meta calculus is a means of modeling mathematical systems of any complexity with building blocks called holons—the holistic units of general systems theory. Mathematically, each holon has a “bandwidth” represented by the dimensionality of its simultaneous solution. In philosophy, the notion of the holon emerged in the 1960s in the writings of Arthur Koestler, and recently the ontology of holons and holarchies has been extensively explored in the integral theory of Ken Wilber. Our thesis is that meta calculus is a mathematical analog of integral theory, and is ontological in the same sense. We use programming examples below to convey the deeper meaning of mathematical holarchies—nested systems of holons.

Meta calculus computational theory emerged during the same period as Koestler’s first publication of the holon, in the software development of rocket propulsion analysis in the Apollo program at NASA and TRW. Nested holonic computation was first used in the technique known variously as state estimation or data assimilation, a computational process of the scientific method. Mathematical models of whole vehicular propulsion systems (CSM and LEM), including the rocket engines, propellant feed systems, and guidance systems, were used to assimilate telemetry data from actual Apollo missions in order to mathematically reconstruct the behavior and performance of these systems for mission analysis. It was the accuracy of this method which validated the payload capability of the Apollo engines as sufficient to carry the Lunar rover vehicle to the moon.

The following table provides a lineage of the languages that evolved from this practice. New languages were necessary to implement the technique of automatic differentiation (AD) arithmetic—a calculus computational architecture that began emerging at TRW and GE, independently during the early 1960s. AD arithmetic automatically computed the partial derivative values of the holon output variables with respect to its input variables, by applying the rules of differential calculus numerically as arithmetic rather than symbolically as algebra. The combination of AD arithmetic and holonic nesting gave rise to the programming paradigm now called meta calculus (AKA synthetic calculus), which completed its evolution in the PROSE language, MC5. PROSE was the first MC language to have the infinitely nestable holonic architecture introduced below. The data assimilation technique had originated in 1958 in assembly-language programs used at TRW for ballistic missile flight reconstruction, and was used to analyze and validate the Mercury and Gemini missions before Apollo. Thus this paradigm evolved over 15 years before it was introduced to the market on supercomputers in 1974.

 

Id

Name

Year

Description

MC1

Model Compiler

1966

NASA AD compiler  (CDC 3600)

MC2

MODTRAN

1967

TRW interpreter/compiler AD (CDC 6600)

MC3

SLANG

1968

TRW Macro Language to MODTRAN (CDC 6600)

MC4

SLANG/CUE

1969

TRW Relocatable Macro Language (CDC 6600)

MC5

PROSE

1972

Batch  Interpretive VM (CDC 6600)

MC6

TSPROSE

1975

Time-shared Interpretive VM (CDC 6400)

MC7

Fortran Calculus

1989

DOS 32-bit protected mode AD compiler

Table A-1  Prior Meta Calculus Language Lineage

Meta calculus resolved and unified higher-mathematical systems engineering into three holonic calculus processes—simulation, correlation, and optimization. These process classes (Si,Co,Op) constitute the “alphabet” of meta calculus (like the A,G,C,T alphabet of DNA). The unified architecture of these mathematical processes is an extension of the subroutine architecture of algebraic programming languages, extending the mathematical scope of holonic subroutines to the inverse math of the scientific method and systems optimization.

Some of the examples given in Chapter 6, and below were originally programmed in SLANG at TRW, and published at the Association of Computing Machinery National Conference in 1969. All of the others were originally programmed in PROSE before 1976. They are presented here in Fortran Calculus (FC), PROSE’s successor, and were executed on a 100MHz 486 processor.

Meta Calculus Holonic Architecture

Scientific-method-level program architecture is composed of dynamic holon structures that are hierarchic execution assemblies associated with the instance of an operator such as FIND, which appears in an upper level operator procedure. This holon operator is associated with the whole operation, and the holon is created when it is executed, according to the parameters appearing in the command. The holon architecture is illustrated in Figure A-1.


Figure A-1 Holon Architecture

 


Conceptually, the holon is a dynamic process of solving a mathematical system subproblem involving simultaneous functions of several variables. The functions are syntactically described in the formulas and subproblem holons contained in the hierarchy of the holon operand—the model procedure. The solution process is controlled by a hidden solver algorithm, specified in the holon operator, which will perform its solution process according to control variables (solver parameters) assigned in a controller procedure, which is also specified in the holon operator if default control variables need be updated. The holon operator calls the solver, and then the solver sets up the problem and executes the model iteratively to solve the holon problem. Before executing the model in each pass, the solver sets triggers in the engine governing the process of automatic differentiation that produces the partial derivative arrays from the model execution. Depending on the solver, partials may only be computed on certain passes. In other passes the model may not be differentiated, or may be differentiated to a different order or even with a different basis vector (coordinate variables).

The holon activation record, called the holog, is the dynamic memory used by the holon. This includes a subset of the coordinate variables (input) and the function variables (output) defined in the holon operator. These holog variables are a subset of the variables involved in the model procedure. Other variables not included in the holog include input parameters which are constants in the holon problem, and auxiliary dependent (output) variables computed in the model but not directly involved in the solution process.

In the following example programs, the holon operator statements are displayed in bold type, and the code examples have been graphically rendered to indicate the deep mathematical sub-structure processes that are been executed “under the hood” in the language semantics. Figure A-0 shows the holon class alphabet symbols. Each holon class is considered a “pantheon” of solver engines, of which only two in each class are mentioned. Solver engines are poly-algorithm systems which operate differently according to the type of mathematical systems being solved. For example, the solver AJAX uses a Newton-Raphson algorithm for determined equation systems, and a Newton-Gauss algorithm to solve  overdetermined equation systems. This is illustrated in the first two examples presented below.


Figure A-0 Holon Graphics Symbols and Associated Engine Examples

In addition, dashed boxes have been added to show the context of automatic differentiation arithmetic, in which all output variables appearing on the left side of formulas are being differentiated with respect to parameters appearing in the FIND statements. In some cases the algorithms within these contexts are being differentiated as well—a feature that is not possible to achieve via formal symbolic differentiation. In cases where the differentiation contexts are nested, then a recursive differential geometry coordinate transformation (not illustrated) is employed to de-nest the differentiation and convert derivatives of the inner nest into derivatives of the outer nest.

Basic Holonic Templates

The following example problems, solved in Fortran Calculus, consist of single level holons, like the AC motor design optimization problem given in Chapter 6. These two problems illustrate two forms of correlation (root finding) of nonlinear equations, which are solved by the FIND … TO MATCH holon operator, where the MATCH phrase identifies equality constraint variables are to be driven to zero within an acceptable tolerance, thus satisfying all of the equations. Both of these problems utilize the Fortran Calculus solver AJAX.

Chemical Equilibrium – Nonlinear Equations

This example is the solution of seven simultaneous nonlinear equations representing a classical chemical equilibrium problem. The model equations simply define a vector of equality constraints to be matched to zero.

 

 

   PROBLEM REACT(10000,1000,1000)

     DIMENSION X(7),F(7)

     DATA X/.5,0.0,0.0,.5,0.0,.5,2.0/

     FIND X; IN MOLES(X,F); BY AJAX; TO MATCH F

   END

 

   MODEL MOLES(X,F)

     DIMENSION X(7),F(7)

     F(1)=X(1)/2+X(2)+X(3)/2-X(6)/X(7)

     F(2)=X(3)+X(4)+2*X(5)-2/X(7)

     F(3)=X(1)+X(2)+X(5)-1/X(7)

     F(4)=-28837*X(1)-139009*X(2)-78213*X(3)+18927*X(4)

  *       +8427*X(5)+13492/X(7)-10690*X(6)/X(7)

     F(5)=X(1)+X(2)+X(3)+X(4)+X(5)-1

     F(6)=400*X(1)*X(4)**3-1.7837E5*X(3)*X(5)

     F(7)=X(1)*X(3)-2.6058*X(2)*X(4)

   END

 


Figure A-1  Chemical Equilibrium Program in Fortran Calculus

The AJAX solver used in this program employs a damped Newton-Raphson algorithm for a determined system of equations (having the same number of equations as unknowns). This problem falls into the meta-calculus category of correlation, in which determined equation systems are considered the “degenerate case” of correlation. Each equation constraint (F vector elements) is matched to within a small tolerance of zero (see Fig. A-2).  Ten Newton iterations were involved in this solution, each computing the Jacobian partial derivative matrix by automatic differentiation arithmetic. In between these iterations, in certain cases, this algorithm executes the model in “damping” sub-iterations of ordinary arithmetic in which partial derivatives are not computed. Note the absence of any output statements in the program, as the solver generates its own report.

 

 

--- AJAX SUMMARY, INVOKED AT REACT[4] FOR MODEL MOLES ----            

 

    CONVERGENCE CONDITION AFTER 10 ITERATIONS

    UNKNOWNS CONVERGED

    CONSTRAINTS SATISFIED 

    ALL SPECIFIED CRITERIA SATISFIED

 

    LOOP NUMBER ....[INITIAL]     9             10

     UNKNOWNS

      X   (  1)  5.000000E-01  3.228708E-01  3.228708E-01

      X   (  2)  0.000000E+00  9.223544E-03  9.223544E-03

      X   (  3)  0.000000E+00  4.601709E-02  4.601709E-02

      X   (  4)  5.000000E-01  6.181717E-01  6.181717E-01

      X   (  5)  0.000000E+00  3.716851E-03  3.716851E-03

      X   (  6)  5.000000E-01  5.767154E-01  5.767154E-01

      X   (  7)  2.000000E+00  2.977863E+00  2.977863E+00

     CONSTRAINTS

      F   (  1)  0.000000E+00 -1.383576E-08 -1.619538E-13

      F   (  2) -5.000000E-01 -9.596910E-08 -1.390105E-11

      F   (  3)  0.000000E+00 -4.798455E-08  1.419820E-11

      F   (  4) -8.815000E+02  4.995033E-04  2.473498E-07

      F   (  5)  0.000000E+00  0.000000E+00  1.079581E-11

      F   (  6)  2.500000E+01  6.865512E-07  1.952332E-08

      F   (  7)  0.000000E+00 -1.060134E-11  9.205562E-12

 

---END OF LOOP SUMMARY

   ELAPSED TIME = 0.61 SECONDS


Figure A-2 Output of Chemical Equilibrium Problem

Admittance Circuit Fitting – Scientific Method

Next is a simple example of the scientific method, or phenomenological model fitting. In this case the model is a set of nonlinear algebraic equations describing the admittance (reciprocal of resistance) derived from Ohm’s law for a particular electric circuit. This problem was originally posed by a PROSE user from Varian Corporation in 1974. It is interesting because it demonstrates the advantage of the scientific method as opposed to ordinary curve fitting. Since the scientific phenomenon is directly represented by the model, the AJAX solver fits the phenomenological discontinuity in the equation as shown in the graph of Figure A-4.  Non-phenomenological equations, such as the polynomials used in ordinary regression would probably not even converge to a solution because of this spike discontinuity in the data.

 

    

      GRAPHICS BOTH

      PROBLEM CIRCUIT (5000,1000,1000)

        COMMON/PARAMS/EL,ELS,CS

        DIMENSION F(21),Y(21),R(21),W(21),YC(21)

        DATA Y/1.273,1.278,1.392,1.604,1.708,1.950,2.148,2.297,

     ~        2.503,2.893,3.305,4.005,5.077,8.069,14.84,39.47,

     ~        -15.06,-9.705,-7.368,-5.286,-3.907/

        EL=1.1E-10 : ELS=-1.1E-10 : CS=1.1E-12

        DO 10 I=1,21

          F(I)=2850+50*I

          W(I)=6.28318E6*F(I)

  10    CONTINUE

        @AXES('CIRFIT',F,Y,YC,W)

        FIND EL,ELS,CS; IN FIT(F,Y,R,W,YC);BY AJAX; TO MATCH R

        @FINISH('CIRFIT',F,YC,W)

      END

 

      MODEL FIT(F,Y,R,W,YC)

        COMMON/PARAMS/EL,ELS,CS

        DIMENSION F(*),Y(*),R(*),W(*),YC(*)

        DO 10 I=1,21

          YC(I)=-1/(W(I)*EL-1/(W(I)*CS-1/(W(I)*ELS)))

          R(I)=1/YC(I)-1/Y(I)

   10   CONTINUE

      END


Figure A-3 Admittance Circuit Program

In this problem, the admittance Y is given in the DATA statement as a set of 21 measurements at 21 frequencies from an unknown electric circuit in which it is desired to duplicate the circuit’s output with a known circuit containing two inductors and a capacitor. The problem is to find the strength of the two inductors and the capacitor to duplicate the unknown circuit.  So, using the admittance formula which is a function of the two unknown inductances (EL and ELS) and the unknown capacitance (CS), we can generate 21 equations corresponding to the 21 measurements in terms of the three unknowns. Thus we have an over-determined set of nonlinear equations. AJAX automatically employs a Newton-Gauss nonlinear least squares approach to solve over-determined systems by minimizing the Euclidean norm of the equations. The constraint vector to be matched in this case is the residual resistance vector R. In this program two graphics subroutines (AXES and FINISH, not shown) are used to plot the admittance points and draw the curve that fits them. The graphical output is shown in Figure A-4


Figure A-4 Admittance Circuit Fit

 

 


--- AJAX SUMMARY, INVOKED AT CIRCUIT[13] FOR MODEL FIT ----           

 

    CONVERGENCE CONDITION AFTER  8 ITERATIONS

      UNKNOWNS CONVERGED

      CONSTRAINTS UNSATISFIED

      ALL SPECIFIED CRITERIA SATISFIED

 

   LOOP NUMBER ...... [INITIAL]       7             8

    UNKNOWNS 

     EL              1.100000E-10  1.497660E-10  1.497660E-10

     ELS            -1.100000E-10 -3.693787E-10 -3.693787E-10

     CS              1.100000E-12  7.472066E-12  7.472069E-12

    OBJECTIVE

     ||G||           2.274699E+00  8.901396E-02  8.901359E-02

 

---END OF LOOP SUMMARY

 

ELAPSED TIME =    8.41 SECONDS

 


Figure A-5 AJAX Summary for Admittance Circuit

 

Nested Holonic Templates – Design Patterns

Nested holons often result in the automatic differentiation of algorithms, a feat which is not feasible with symbolic differentiation, and the reason I have always considered tools like MACSYMA as limited in regard to numerical computing. They have great usefulness in higher-order model transformation, and I have high hopes for their use in automating the modeling process, particularly in regard to modeling systems of partial-differential equations (PDEs).

Pilot Ejection Altitude-Speed Profiling

During the mainframe era of problem-oriented languages before 1980, the pilot-ejection problem was a favorite model used to benchmark the syntax of continuous-system simulation languages used for solving simultaneous ordinary differential equations (ODEs), such as the CSSL standard, from which I took the equations of this problem. I programmed them in PROSE as two different examples in its manual. The first was a simulation approach, like the CSSL version, which computed the altitude-speed profile of safe ejection by trial and error as an interpolation study of many simulation runs.  A user would have to run each trajectory individually in this case.

The second one, which this FC example (Fig. A-6) duplicates, was a pedagogical demonstration of the power of automatic differentiation in PROSE by nesting the ODE simulation holon inside an AJAX holon to solve the problem as a repeated two-point boundary-value problem computing and plotting the altitude-speed safe-ejection profile in a single run (Fig A-8). Six two-point boundary-value problems are solved for a range of altitudes (0-50,000 ft) by placing the FIND statement for the AJAX holon in a loop. The aircraft speed and seat trajectory time are solved for as unknowns by the FIND statement to match the boundary condition constraints GX and GY which specify minimum clearance for safe ejection.

The differential equations reside in the model MOTION. The initial-value problem to solve these equations is defined and initialized by the INITIATE statement, which creates the ISIS solver holon. It identifies the rate variables to be integrated (THEDOT, VDOT, XDOT, and YDOT), the corresponding state variables that are the output of integration (THETA, V, X, and Y),  the independent variable T,  the integration step size DT, and the upper limit of integration TP.

 

 

      GLOBAL ALL

      GRAPHICS FILE

      PROBLEM EJECT  ! Pilot Ejection Profile

        DIMENSION ALT(6),VEL(6)

        CHARACTER FCSINT*2,PLANE*7

        DATA ALT/0,10000,20000,30000,40000,50000/

        SMASS=7 : G=32.2 : CD=1 : VE=40 : THETAD=15 : TIME=1.0

        G=10 : Y1=4 : VA=100

        DO 10 I=1,6

          H=ALT(I) : IT=0 : PLANE='PLANE'//FCSINT(I)

          @AXES(PLANE,'PILOT EJECTION')

          IF(H.LE.35332) THEN

             RHO=0.002378*(1-.689E-5*H)**4.256

          ELSE

             RHO=0.00315/EXP(1.452+(H-35332)/20950)

          ENDIF

         FIND VA,TIME; IN SEAT(PLANE);

      &     BY AJAX(ACON); TO MATCH GX,GY

CORRELATION          @DISPLAY(PLANE)

          VEL(I)=VA

   10   CONTINUE

        @SAFE('PROFILE')

      END

Automatic Differentiaton      MODEL SEAT(PLANE)

        CHARACTER*2 FCSINT,NI,PLANE*7

        VX=SQRT(VA**2)-VE*SIND(THETAD) : VY=VE*COSD(THETAD)

        V=SQRT(VX*VX+VY*VY) : THETA=ATAN(VY/VX)

        X=0 : Y=Y1 : T=0 : DT=ABS(TIME)/20 : DP=4*DT : TP=T+DP

        IT=IT+1 : NI=FCSINT(IT)

        IF(IT.GT.9) THEN

          @POINT(PLANE,'I'//NI,X,Y)

        ELSE

          @CURVE(PLANE,'I'//NI,X,Y)

        ENDIF

       INITIATE ISIS; FOR MOTION; EQUATIONS

     &   THEDOT/THETA, VDOT/V, XDOT/X, YDOT/Y; OF T; STEP DT; TO TP

        DO WHILE (T.LT.TIME)

         INTEGRATE MOTION; BY ISIS

          @CURVE(PLANE,'I'//NI,X,Y)

          TP=TP+DP

SIMULATION        END DO

        GX=-X-30     ! Boundary condition on X at T=TIME

        GY=Y-20      ! Boundary condition on Y at T=TIME

       TERMINATE MOTION

      END

      MODEL MOTION  ! Differential Equations

        D=0.5*RHO*CD*S*V*V  ! Drag equation

        THEDOT=-G*COS(THETA)/V

        VDOT=-D/SMASS-G*SIN(THETA)

        XDOT=V*COS(THETA)-ABS(VA)

        YDOT=V*SIN(THETA)

      END

 


Figure A-6  Pilot Ejection Safe Profile Program

 

The INTEGRATE statement then generates the trajectory using the solver ISIS. In this case the integration interval is a segment of the trajectory between plotted points as defined by the upper limit TP which is advanced forward by four integration steps in each loop. In this example, the numerical integration solver algorithm ISIS of the inner holon is automatically differentiated as part of the model of the outer holon. Thus each seat-trajectory curve in the demo graphs of Fig A-7 represents an iteration of AJAX in the solution of the two-point boundary-value problem, which finally converges to the trajectory which hits the target (the minimum safe trajectory). Integration step size DT is dependent upon the TIME parameter computed by AJAX. This dependence is the key linkage for propagation of partial derivatives which enables the nested solution to converge.


 Figure A-7(a) Pilot Ejection Demo Output for Zero Altitude 

 



Figure A-7(b) Pilot Ejection Demo Output for 10,000 feet

 


   Figure A-7(c) Pilot Ejection Demo Output for 20,000 feet


Figure A-7(d) Pilot Ejection Demo Output for 30,000 feet


    Figure A-7(e) Pilot Ejection Demo Output for 40,000 feet


     Figure A-7(f) Pilot Ejection Demo Output for 50,000 feet


Figure A-8 Safe Ejection Altitude-Speed Profile

When this pedagogical demo problem is performed interactively, the program pauses for each of the Figure A-7 plots to give the user a visual sense of the intermediate solutions of the equations. In the original CSSL version of this model it was necessary to do this for each seat trajectory—involving many runs. In this boundary value version, however, the total solution to the whole problem (Figure A-8) could be computed in a single run. In a non-pedagogical scenario the program could be executed as a batch process, generating all the graphs as output. The output reports generated by AJAX are not shown in this case.

Bipropellant Rocket – Implicit Differential Equations

This example (Fig. A-9) continues the discussion of Implicit Differential Equations in Chapter 6 (See Fig. 6-14). This involves nesting a correlation (root-finding) holon inside a simulation holon. This model is a simplification of the model of the oxidizer and fuel flow rates of the Apollo CSM service-module rocket engine, where a longitudinal (pogo) vibration influence is added to the flow-rate equations. The inner correlation holon solves the implicit flow rate equations during each integration step of the outer simulation holon. The simulation solver, JANUS, an Adams-Moulton predictor-corrector algorithm, calls the model IDEQ multiple times during each integration step to sample the rate variables XDOT and YDOT solved by the Newton–Raphson algorithm of the solver AJAX. Fig. A-10 shows the propellant weight and flow-rate curves (scaled and biased to show all four curves on the same graph) from the solution.

 

      GLOBAL ALL    

      PROBLEM IMPDES

        X=14000 : Y=7000       ! Initial Conditions

        XDOT=-50 : YDOT=-25    ! Initial Rate Guesses

        T=0 : DT=.5 : TP=2 : TF=TP

        @AXES(‘PLOT’)

        INITIATE JANUS; FOR IDEQ EQUATIONS

     ~     XDOT/X, YDOT/Y; OF T; STEP DT; TO TF

        DO WHILE (TF.LE.50)

          INTEGRATE IDEQ; BY JANUS

          @CURVES(‘PLOT’)

          TF=TF+TP

SIMULATION        END DO

      END

 

      MODEL IDEQ   !    Implicit Differential Equations

        FIND XDOT, YDOT; IN IRATE; BY AJAX(ACON); TO MATCH GX,GY

      END

 

      CONTROLLER ACON(AJAX)

        SUMMARY=0

CORRELATION      END

 


Automatic Differentiaton      MODEL IRATE   !   Implicit Rate Equations

        GX=XDOT+3.2*SQRT(1-(XDOT+YDOT)*

     ~     (1.15+57.5/(20000+X+Y)))*(1+.1*EXP(-T/10)*SIN(1.5708*T))

        GY=YDOT+1.6*SQRT(1-(XDOT+YDOT)*

     ~     (1.15+36.2/(20000+X+Y))*(1+.1*EXP(-T/10)*SIN(1.5708*T))

      END

 


Figure A-9  Implicit Differential Equations Program


Figure A-10 Implicit Differential Equations Output

 


Three-Stage Rocket Design Optimization

This is an example of nesting a correlation (AJAX) holon inside an unconstrained (HERA) optimization holon, to solve an equality-constrained optimization problem of a three-stage rocket. This approach to equality-constrained optimization forces the search process to follow the constraint surface in finding the extremum point of the optimization. This problem was originally solved in 1966 using the Restricted Optimization Program (ROP), a second-order Newton optimizer which became the API of MODTRAN (MC2). The model was input to the ROP program using what might be considered MC0—the  “assembly” form of AD arithmetic. The rocket model was a preliminary design of a Mars surface sample and retrieval vehicle in an early TRW study for JPL.

 

 

      GLOBAL ALL

      PROBLEM ROCKET(20000,2000,2000)

        DIMENSION SPI(3),SPIVAC(3),TBURN(3),THRUST(3),XIP(3),WPROP(3)

        DIMENSION  RATIO(3),WSTAGE(3),STRFAC(3),DELV(3),G(2)

        DATA THRUST /350,1500,4100/,TBURN/110,100,180/

        DATA  XIP/5E-3,0,0/,SPIVAC/3*315/

        FIND THRUST(1),THRUST(2),TBURN(2),TBURN(3); IN STAGES;

     &     BY HERA;  REPORTING DLVTOT,TBTOT;

     &     TO MINIMIZE WEIGHT    

      END

 

      MODEL STAGE

        FIND THRUST(3), TBURN(1); IN EQS; BY AJAX; TO MATCH G

      END

 

      MODEL EQS

        ! THREE STAGE ROCKET PHENOMENA MODEL   

        DATA GC,WPAYLD,DELVIP,TBIP/32.174,50,2.8E4,400/

        DLVTOT=0 : TBTOT=0

        WEIGHT=WPAYLD

        DO I=1,3

          SPI(I)=SPIVAC(I)*(1-XIP(I))

          WPROP(I)=THRUST(I)*TBURN(I)/SPI(I)

          WSTAGE(I)=0.0234*THRUST(I)+WPROP(I)+1.255*WPROP(I)**0.704+4

          WEIGHT=WEIGHT+WSTAGE(I)

          RATIO(I)=WEIGHT/(WEIGHT-WPROP(I))

          DELV(I)=GC*SPI(I)*LOG(RATIO(I))

          DLVTOT=DLVTOT+DELV(I)

          TBTOT=TBTOT+TBURN(I)

        END DO

        G(1)=DLVTOT-DELVIP ! Total delta V constraint

        G(2)=TBTOT-TBIP    ! Total burn time constraint

      END

 


      Figure A-11 Three Stage Rocket Design Optimization

Of particular interest in this type of nested holonic structure is that both holons are implicit, thus the propagation of partial derivatives from  the inner holon to the outer holon is not simply the differentiation of an algorithm, as in the case of the pilot-ejection problem shown above. In this case, upon convergence of the inner holon, a differential geometry coordinate transformation is required to convert the inner holon derivatives into the coordinate system of the outer holon. This is a recursive procedure that is inherent in the architecture of the meta-calculus engine. It is a wrapper harness for all of the implicit solvers that propagate partial derivatives.

 

 

--- AJAX SUMMARY, INVOKED AT STAGES[19] FOR MODEL EQNS ----           

   CONVERGENCE CONDITION AFTER  3 ITERATIONS

      UNKNOWNS CONVERGED

      CONSTRAINTS SATISFIED

      ALL SPECIFIED CRITERIA SATISFIED

 LOOP NUMBER ......   [INITIAL]      1             2             3

 UNKNOWNS

   THRUST  (  3)    4.100000E+03  4.435707E+03  4.454498E+03  4.454535E+03

   TBURN   (  1)    1.100000E+02  1.200000E+02  1.200000E+02  1.200000E+02

 CONSTRAINTS

   G       (  1)   -4.907225E+02 -2.130769E+01 -4.224610E-02 -1.666122E-07

   G       (  2)   -1.000000E+01  0.000000E+00  0.000000E+00  0.000000E+00

---END OF LOOP SUMMARY

--- AJAX SUMMARY, INVOKED AT STAGES[19] FOR MODEL EQNS ----           

   CONVERGENCE CONDITION AFTER  3 ITERATIONS

      UNKNOWNS CONVERGED

      CONSTRAINTS UNSATISFIED

      ALL SPECIFIED CRITERIA SATISFIED

 LOOP NUMBER ......   [INITIAL]      1             2             3

 UNKNOWNS

   THRUST  (  3)    4.454535E+03  5.067309E+03  5.107165E+03  5.107312E+03

   TBURN   (  1)    1.200000E+02  1.206324E+02  1.206324E+02  1.206324E+02

 CONSTRAINTS

   G       (  1)   -6.753633E+02 -3.875873E+01 -1.420019E-01 -1.917695E-06

   G       (  2)   -6.324392E-01  0.000000E+00  0.000000E+00  0.000000E+00

 LOOP NUMBER .........   [INITIAL]         3

---END OF LOOP SUMMARY

--- AJAX SUMMARY, INVOKED AT STAGES[19] FOR MODEL EQNS ----           

   CONVERGENCE CONDITION AFTER  2 ITERATIONS

      UNKNOWNS CONVERGED

      CONSTRAINTS UNSATISFIED

      ALL SPECIFIED CRITERIA SATISFIED

 LOOP NUMBER ......   [INITIAL]      1             2

 UNKNOWNS

   THRUST  (  3)    5.107312E+03  5.102872E+03  5.104214E+03

   TBURN   (  1)    1.206324E+02  1.159654E+02  1.159654E+02

 CONSTRAINTS

   G       (  1)    3.696945E+01 -1.292751E+00 -1.572368E-04

   G       (  2)    4.667031E+00  0.000000E+00  0.000000E+00

---END OF LOOP SUMMARY

--- AJAX SUMMARY, INVOKED AT STAGES[19] FOR MODEL EQNS ----           

   CONVERGENCE CONDITION AFTER  2 ITERATIONS

      UNKNOWNS CONVERGED

      CONSTRAINTS SATISFIED

      ALL SPECIFIED CRITERIA SATISFIED

 LOOP NUMBER ......   [INITIAL]      1             2

 UNKNOWNS

   THRUST  (  3)    5.104214E+03  5.110939E+03  5.110944E+03

   TBURN   (  1)    1.159654E+02  1.158542E+02  1.158542E+02

 CONSTRAINTS

   G       (  1)   -5.639499E+00 -4.708355E-03 -2.080924E-09

   G       (  2)    1.112060E-01  0.000000E+00  0.000000E+00

---END OF LOOP SUMMARY

---- HERA SUMMARY, INVOKED AT ROCKET[16] FOR MODEL STAGES ----        

   CONVERGENCE CONDITION AFTER  3 ITERATIONS

      UNKNOWNS NOT CONVERGED

      OBJECTIVE CRITERION SATISFIED

      ALL SPECIFIED CRITERIA SATISFIED

 LOOP NUMBER ......   [INITIAL]      1             2             3

 UNKNOWNS

   THRUST  (  1)    3.500000E+02  3.135112E+02  3.271673E+02  3.282804E+02

   THRUST  (  2)    1.500000E+03  1.279613E+03  1.278271E+03  1.277462E+03

   TBURN   (  2)    1.000000E+02  1.262872E+02  1.326462E+02  1.330776E+02

   TBURN   (  3)    1.800000E+02  1.530804E+02  1.513884E+02  1.510682E+02

 OBJECTIVE

   WEIGHT           3.814512E+03  3.780665E+03  3.778655E+03  3.778647E+03

 OTHER VARIABLES

   DLVTOT           2.800000E+04  2.800000E+04  2.800000E+04  2.800000E+04

   TBTOT            4.000000E+02  4.000000E+02  4.000000E+02  4.000000E+02

---END OF LOOP SUMMARY

ELAPSED TIME =    0.11 SECONDS

 


        Figure A-12 Three Stage Rocket Design Optimization Output

 

Chemical Kinetics – Scientific Method

This is another scientific method template with a set of differential equations as model, in which sampled data at periodic intervals over time is to be fitted by the output of each equation. It is therefore a multi-point boundary-value problem which solves for the parameters of the differential equations required to fit the data. Since the differential equations in this problem are chemical reaction rate equations, this is a classical problem in chemical kinetics, where it is dealing with non-equilibrium chemical transients—unlike the chemical equilibrium problem above (idealized reactions with infinite reaction rates).

The Apollo spacecraft rocket engines had to deal with non-equilibrium chemical kinetics, whereas the Titan rockets, which used the same propellants, expanded to higher ranges of pressure, where equilibrium chemistry models were accurate predictors of rocket performance.  Models such as this, having many differential equations were used to estimate and verify performance forecasts from propellant testing during the development of the Apollo engines.

This same design pattern template is well known in the pharmaceutical industry, where it is called “pharmaco-kinesis”—determining the parameters of the reaction rate equations by periodic measurement of the concentrations of reaction constituents in the blood, and fitting these measurements to the output of integrating the differential equations.

The example problem consists for four differential equations for the varying concentrations of four chemical constituents, A, B, C, and D. Samples of A and C concentrations spread over the time interval are the data to be matched. Thus a simulation holon in model CURFIT (Fig A-13) integrates the four differential equations in model REACTION. This ATHENA holon is nested inside an unconstrained optimization holon using the solver HERA which minimizes the sum-squared error. In this case the ATHENA holon is established once at the top level of the program (outside the optimization loop) by the INITIATE statement preceding the FIND statement.

For pedagogical illustration purposes like the pilot-ejection problem, this program generates a graph during each optimization iteration, showing the constituent curves and the corresponding measurement data, resulting in the six graphs of Fig. A-14. The demo version thus pauses when each graph is generated so the user can gain a sense of the intermediate trajectory solutions of the differential equations, as the optimizer HERA is computing the intermediate values of the parameters which govern the trajectories. In a non-pedagogical “batch” version, only the final graph, Figure A-14(g), and output report, Figure A-15, would be generated.

 

 

      GLOBAL ALL

      GRAPHICS FILE

      PROBLEM CHEMPARE  ! Chemical Kinetics Parameter Estimation

        DIMENSION TM(6),CM(6),AM(6)

        DATA TM/10,20,30,40,50,60/  ! Time points for measurements

        DATA CM/0.419,0.563,0.629,0.666,0.689,0.708/  ! C samples

        DATA AM/0.483,0.281,0.191,0.134,0.097,0.065/  ! A samples

        DATA NM,A0,B0,C0,D0,DT/6,1.0,1.03,0.0,0.0,2.0/! Initial Cond.

        DATA B1,B2/0.015,0.015/   ! Optimization step bounds

        P1=0.01: P2=0.05          ! Parameter starting guesses

       INITIATE ATHENA; FOR REACTION;

     &    EQUATIONS DADT/A,DBDT/B,DCDT/C,DDDT/D;

     &    OF T; STEP DT; TO TF

       FIND P1,P2; IN CURFIT;

     &    BY HERA(SET); WITH BOUNDS B1,B2; TO MINIMIZE ERROR

      END

 

      MODEL CURFIT

Automatic Differentiaton        CHARACTER*2 FCSINT,NI

        DATA IT/0/ : NI=FCSINT(IT)

        @AXES('ITER'//NI,'ITERATION '//NI)

        A=A0: B=B0: C=C0: D=D0: T=0: TF=0: ERROR=0: I=1

        DO WHILE (I.LE.NM)

          TF=TF+DT

         INTEGRATE REACTION; BY ATHENA

          @CURVES('ITER'//NI)

          IF(TF.EQ.TM(I)) THEN

            ERROR=ERROR+(AM(I)-A)**2+(CM(I)-C)**2  ! Cumulative error

            @MEASURED('ITER'//NI,AM(I),CM(I),TM(I)) ! Plot measurement

            I=I+1

          ENDIF

        END DO

        @SHOW('ITER'//NI)       ! Show graph for this iteration

        IT=IT+1

      END

 

      MODEL REACTION    ! Rates of reaction differential equations

        DCDT=P1*A*B

        DADT=-(DCDT+(.01+P2)*A

        DBDT=-(DCDT+.05*B*D)

        DDDT=DBDT-DADT

      END

 

      CONTROLLER SET(HERA)

         DELTA=1E-2:DETAIL=1: ADJUST=2

      END

 


Figure A-13 Chemical Kinetics Parameter Estimation Program

 


 



 


 


 


 


 


 


---- HERA SUMMARY, INVOKED AT CHEMPARE[12] FOR MODEL CURFIT ----      

    CONVERGENCE CONDITION AFTER  6 ITERATIONS

      UNKNOWNS CONVERGED

      OBJECTIVE CRITERION UNSATISFIED

      ALL SPECIFIED CRITERIA SATISFIED

  LOOP NUMBER ....[INITIAL]       5             6

  UNKNOWNS

       P1        1.000000E-02  7.227140E-02  7.315935E-02

       P2        5.000000E-02  3.450409E-03  3.565576E-03

  OBJECTIVE

       ERROR     1.597781E+00  1.967029E-04  1.275779E-04

 ---END OF LOOP SUMMARY

 


 Figure A-15  Summary Report for HERA Optimizer

Optimal Design and Control (ODC) Template

This example continues the discussion of the optimal design and control problem template from Chapter 6.  The resulting design-pattern template involves three nested holons, represented by the bold command statements. The inner two comprises the two-point boundary-value problem, using an ISIS simulation holon nested within an AJAX correlation holon, and this in turn is nested within a HERA optimization holon. (Fig. A-16).

 

      GLOBAL ALL

      PROBLEM OPTDES   ! Optimal design & control

        FIND A; IN TWOPT; BY HERA; TO MINIMIZE OBJ

      END

 

      MODEL TWOPT ! Two point boundary value problem model

        Y0=.7 ! Guess for initial condition

        FIND Y0; IN TRAJ; BY AJAX; TO MATCH Y1

        OBJ=Z/2+A**2/2

      END

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 



Figure A-16  Optimal Design and Control Program

Figure A-17  shows the output of this application. It was originally programmed in SLANG, and published in the proceedings of the ACM National Conference in 1969, referencing the origin of the problem from an advanced mathematical text by Richard Bellman and Robert Kalaba, which contained the original FORTRAN program containing an algorithm they invented, called quasilinearization. The FORTRAN program (Fig. 6-13) also used built-in integration solvers not illustrated in the textbook. But the statement of the algorithm in the program, which required 125 statements, made the problem being solved incomprehensible.  The Fortran Calculus version in Figure A-16 illustrates the metaphoric (non-algorithmic) statement of the problem, quite simple even in a procedural language.

 

 

---AJAX SUMMARY, INVOKED AT TWOPT[8] FOR MODEL TRAJ ----

     CONVERGENCE CONDITION AFTER  1 ITERATIONS

       UNKNOWNS CONVERGED

       CONSTRAINTS SATISFIED

       ALL SPECIFIED CRITERIA SATISFIED

   LOOP NUMBER ......  [INITIAL]      1

   UNKNOWNS

       Y0         0.700000E+00 -0.761580E+00

   CONSTRAINTS

       Y1         0.225529E+01 -0.183664E-15

 ---END OF LOOP SUMMARY

 ---AJAX SUMMARY, INVOKED AT TWOPT[8] FOR MODEL TRAJ ----

     CONVERGENCE CONDITION AFTER  1 ITERATIONS

       UNKNOWNS CONVERGED

       CONSTRAINTS SATISFIED

       ALL SPECIFIED CRITERIA SATISFIED

   LOOP NUMBER ......  [INITIAL]      1

   UNKNOWNS

       Y0         0.700000E+00 -0.642454E+00

   CONSTRAINTS

       Y1         0.247563E+01  0.602491E-15

 ---END OF LOOP SUMMARY

 ---AJAX SUMMARY, INVOKED AT TWOPT[8] FOR MODEL TRAJ ----

     CONVERGENCE CONDITION AFTER  1 ITERATIONS

       UNKNOWNS CONVERGED

       CONSTRAINTS SATISFIED

       ALL SPECIFIED CRITERIA SATISFIED

   LOOP NUMBER ...... [INITIAL]      1

   UNKNOWNS

       Y0        0.700000E+00 -0.640252E+00

   CONSTRAINTS

       Y1        0.248091E+01 -0.174665E-15

 ---END OF LOOP SUMMARY

 ---AJAX SUMMARY, INVOKED AT TWOPT[8] FOR   GLOBAL ALL                     

     CONVERGENCE CONDITION AFTER  1 ITERATIONS

       UNKNOWNS CONVERGED

       CONSTRAINTS SATISFIED

       ALL SPECIFIED CRITERIA SATISFIED

   LOOP NUMBER ...... [INITIAL]      1

   UNKNOWNS

       Y0        0.700000E+00 -0.640252E+00

   CONSTRAINTS

       Y1        0.248091E+01 -0.174665E-15

 ---END OF LOOP SUMMARY

 ---AJAX SUMMARY, INVOKED AT TWOPT[8] FOR MODEL TRAJ ----

     CONVERGENCE CONDITION AFTER  1 ITERATIONS

       UNKNOWNS CONVERGED

       CONSTRAINTS SATISFIED

       ALL SPECIFIED CRITERIA SATISFIED

   LOOP NUMBER ...... [INITIAL]      1

   UNKNOWNS

       Y0        0.700000E+00 -0.640251E+00

   CONSTRAINTS

       Y1        0.248091E+01 -0.234730E-15

 ---END OF LOOP SUMMARY

 ---HERA SUMMARY, INVOKED AT OPTDES[3] FOR MODEL TWOPT ----

      CONVERGENCE CONDITION AFTER  3 ITERATIONS

       UNKNOWNS CONVERGED

       OBJECTIVE CRITERION SATISFIED

       ALL SPECIFIED CRITERIA SATISFIED

   LOOP NUMBER ......[INITIAL]      1             2             3

   UNKNOWNS

       A        0.000000E+00  0.228183E+00  0.232900E+00  0.232902E+00

   OBJECTIVE

       OBJ      0.380810E+00  0.347279E+00  0.347265E+00  0.347265E+00

---END OF LOOP SUMMARY

 

  ELAPSED TIME =    5.27 SECONDS

 


  Figure A-17  Summary Output of Optimal Design & Control Problem

 

ODC Template Expanded – Wing Design Optimization

Figure A-18 shows a application which is an expanded form of the ODC template, a ten parameter constrained lift maximization of a cantilever beam differential equation two-point boundary-value problem. The obvious point being made here is that self-organizing mathematical design-pattern templates arise in meta calculus, as they do in object-oriented programming. And because the holon alphabet only has three elements, the variety of design patterns is quite small.

 

 

! Wing Design Optimization

   PROBLEM WING (40000,5000,2000)

     DYNAMIC  XS,Y,ALPHA,ALPHAL,ALPHAU,DELALP,EIS,EIL,EIU,DELEI

     @NODIAG(1006): CALL SETUP

     FIND EIS,ALPHA; ! Flexural rigidity vector, length fraction vector

  ~    IN FLEX; BY THOR(TCON); WITH BOUNDS DELEI,DELALP;

  ~    WITH LOWERS EIL,ALPHAL;  WITH UPPERS EIU,ALPHAU;

  ~    HOLDING WEIGHT; MATCHING ALTOTL;

  ~    TO MAXIMIZE WINGLIFT  ! Lift/Weight

     PRINT *,' SOLUTION TO ODE-S'

     TABULATE  XS,Y

   END

 

   MODEL FLEX

     FIND DY0; IN BEAMIVP;  ! Beam equation boundary value problem

  ~    BY AJAX(KNOB); TO MATCH YPA

     EIAVG=0

     DO K=1,NEI

       EIAVG=EIAVG+ALPHA(K)*EIS(K) ! Avg. flexural rigidity

     END DO

     LFLG=1: WINGLIFT=$INTEGSI(FORCE,O,WL,1D-4)/EIAVG

     LFLG=0: WEIGHT=WTMAX-EIAVG       ! Weight (inequality) constraint

     ALTOTL=ARRAYSUM(ALPHA)-1 ! Total length (equality) constraint

     ROWPRINT  EIAVG,WINGLIFT,WEIGHT,ALTOTL

     TABULATE  EIS,ALPHA

   END    

 

   MODEL BEAMIVP             ! Beam equation initial value problem

     X=O: Y0=0: DY0DX=DY0: XF=H:  Y(1)=Y0

     INITIATE ATHENA(CNTRL); FOR BEAMODE; ! Integrate beam equation

  ~    EQUATIONS D2Y0DX2/DY0DX,DY0DX/Y0; OF X; STEP H; TO XF

     DO I=2,NPTS

       INTEGRATE BEAMODE; BY ATHENA

        Y(I)=Y0:  XF=XF+H

     END DO

     TERMINATE  BEAMODE

     YPA=DY0DX-DYLDX  ! Boundary constraint - must be zero

   END

 

   MODEL BEAMODE  ! Cantilever beam differential equation

     SUM=0

     DO J=1,NEI

       IF(X.GE.SUM*WL.AND.X.LT.(ALPHA(J)+SUM)*WL) THEN

         EI=EIS(J):  GOTO 20

       ENDIF

       SUM=SUM+ALPHA(J)

     END DO

20   D2Y0DX2=-$INTEGSI(FORCE,O,X,4)/EI*(1+DY0DX**2)**1.5

   END

 

   FMODEL FORCE(R) ! Quadrature integrand - lift or moment

     IF(LFLG.GT.0) THEN  ! Calculate lift integrand

       FORCE=A*COS(PI*R/(1.95*WL))-EPSLN*TERP(R)

     ELSE  ! Calculate moment integrand

       FORCE=R*(A*COS(PI*R/(1.95*WL))-EPSLN*TERP(R))

     ENDIF

   END

 

   FMODEL TERP(Z) ! Lagrange interpolation of wing formula

     IF(Z.EQ.XS(1)) THEN

        TERP=Y(1)

     ELSEIF(Z.EQ.XS(NPTS)) THEN

        TERP=Y(NPTS)

     ELSEIF(Z.LT.XS(2)) THEN

        TERP=POLY(XS(1),XS(2),XS(3),Y(1),Y(2),Y(3),Z)

     ELSE

       DO II=2,NPTS-1

         IF(Z.GE.XS(II)) GOTO 20

       END DO 

20     TERP=POLY(XS(II-1),XS(II),XS(II+1),Y(II-1),Y(II),Y(II+1),Z)

     ENDIF

   END

 

   FMODEL POLY(XI,XJ,XK,YI,YJ,YK,ZZ)

     GI=(ZZ-XJ)*(ZZ-XK)/((XI-XJ)*(XI-XK))

     GJ=(ZZ-XI)*(ZZ-XK)/((XJ-XI)*(XJ-XK)): GK=(ZZ-XI)*(ZZ-XJ)/((XK-XI)*(XK-XJ))

     POLY=GI*YI+GJ*YJ+GK*YK

   END

 

 

    PROCEDURE SETUP

     PI=3.141592654:O=0: WL=1: A=6: EPSLN=0.1: DY0=0.01: DYLDX=0:

     H=0.05: NEI=5: WTMAX=7:NPTS=(WL-O)/H+1

     ALLOT ALPHA(NEI),ALPHAL(NEI),ALPHAU(NEI),ALPHAU(NEI),DELALP(NEI)

     ALLOT EIS(NEI),EIL(NEI),EIU(NEI),DELEI(NEI): 

     ALLOT XS(NPTS),Y(NPTS)

     <<ALPHA>>=(0.2): <<ALPHAL>>=(0.1): << ALPHAU>>=(1.0)

     <<DELALP>>=(0.0333): <<DELEI>><EI>=(0.333): <<EIS>>=DATA(6,5,4,3,2)

     <<EIL>>=DATA(3,2,2.5,1,1):  <<EIU>>=DATA(10,7,5,4,4)

     DO I=1,NPTS

       XS(I)=O+(I-1)*H:  Y(I)=EPSLN*SIN(PI*XS(I)/2)

     END DO       

     TABULATE EIL,EIU,ALPHA<EIU>L

   END

  

   CONTROLLER TCON(THOR)

     REMAX=12

     DETAIL=1

   END

 

   CONTROLLER KNOB(AJAX)

     REMAX=10

     DAMP=0

     ZERO=H**3

   END

 

   CONTROLLER CNTRL(ATHENA)

     ORDER=2

   END


Figure A-18(a)  Wing Design Optimization Holarchy

 

Figure A-19 provides the output of the wing design solution process.

 

---THOR SUMMARY, INVOKED AT WING[6] FOR MODEL FLEX ----               

    CONVERGENCE CONDITION AFTER  5 ITERATIONS

      UNKNOWNS CONVERGED

      OBJECTIVE CRITERION SATISFIED

      ALL SPECIFIED CRITERIA SATISFIED

 LOOP NUMBER ......[INITIAL]       1             2             3

 UNKNOWNS

   EIS    (  1)   6.000000E+00  5.334000E+00  4.668000E+00  4.002000E+00

   EIS    (  2)   5.000000E+00  4.334000E+00  3.668000E+00  3.002000E+00

   EIS    (  3)   4.000000E+00  3.334000E+00  2.668000E+00  2.500000E+00

   EIS    (  4)   3.000000E+00  2.334000E+00  1.668000E+00  1.002000E+00

   EIS    (  5)   2.000000E+00  1.334000E+00  1.000000E+00  1.000000E+00

   ALPHA  (  1)   2.000000E-01  1.750000E-01  1.084000E-01  1.000000E-01

   ALPHA  (  2)   2.000000E-01  1.700000E-01  1.034000E-01  1.000000E-01

   ALPHA  (  3)   2.000000E-01  1.625000E-01  1.625000E-01  1.000000E-01

   ALPHA  (  4)   2.000000E-01  2.259000E-01  2.925000E-01  3.002000E-01

   ALPHA  (  5)   2.000000E-01  2.666000E-01  3.332000E-01  3.998000E-01

 OBJECTIVE

   WINGLIFT       9.278165E-01  1.197616E+00  1.729176E+00  2.239101E+00

 INEQUALITY CONSTRAINTS

   WEIGHT         3.000000E+00  3.905100E+00  4.860078E+00  5.349000E+00

 EQUALITY CONSTRAINTS

   ALTOTL         0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00

 

 LOOP NUMBER ......[INITIAL]       4             5

 UNKNOWNS

   EIS    (  1)   6.000000E+00  3.336000E+00  3.000000E+00

   EIS    (  2)   5.000000E+00  2.336000E+00  2.000000E+00

   EIS    (  3)   4.000000E+00  2.500000E+00  2.500000E+00

   EIS    (  4)   3.000000E+00  1.000000E+00  1.000000E+00

   EIS    (  5)   2.000000E+00  1.000000E+00  1.000000E+00

   ALPHA  (  1)   2.000000E-01  1.000000E-01  1.000000E-01

   ALPHA  (  2)   2.000000E-01  1.000000E-01  1.000000E-01

   ALPHA  (  3)   2.000000E-01  1.000000E-01  1.000000E-01

   ALPHA  (  4)   2.000000E-01  2.336000E-01  2.336000E-01

   ALPHA  (  5)   2.000000E-01  4.664000E-01  4.664000E-01

 OBJECTIVE

   WINGLIFT       9.278165E-01  2.436545E+00  2.549457E+00

 INEQUALITY CONSTRAINTS

   WEIGHT         3.000000E+00  5.482800E+00  5.550000E+00

 EQUALITY CONSTRAINTS

   ALTOTL         0.000000E+00  0.000000E+00  0.000000E+00

 

 SOLUTION TO ODE-S

                  XS           Y

        1.  0.000000E+00  0.000000E+00

        2.  5.000000E-02  2.452446E-02

        3.  1.000000E-01  4.904389E-02

        4.  1.500000E-01  7.353753E-02

        5.  2.000000E-01  9.797026E-02

        6.  2.500000E-01  1.223121E-01

        7.  3.000000E-01  1.465227E-01

        8.  3.500000E-01  1.703469E-01

        9.  4.000000E-01  1.935454E-01

       10.  4.500000E-01  2.159355E-01

       11.  5.000000E-01  2.373134E-01

       12.  5.500000E-01  2.574549E-01

       13.  6.000000E-01  2.761171E-01

       14.  6.500000E-01  2.930404E-01

       15.  7.000000E-01  3.079527E-01

       16.  7.500000E-01  3.205759E-01

       17   8.000000E-01  3.306368E-01

       18   8.500000E-01  3.378819E-01

       19.  9.000000E-01  3.420975E-01

       20.  9.500000E-01  3.431316E-01

 

 ELAPSED TIME =   32.35 SECONDS

 


   Figure A-19 Output of Wing Design Optimization Program

Systems Engineering End-User Perspective

These examples demonstrate a programming perspective that is directly focused on application programming in science and engineering,  where the end-user is a casual programmer who uses computers to solve sophisticated scientific problems, but is not necessarily a software expert, and doesn’t have time to become one. Although the programs shown in this Appendix are textbook-size problems, that does not imply any limitations on this technology.  It is only limited by the CPU power and memory capacity of computers, which has magnified by many orders of magnitude since these programs were originally written. The computer in which they were run (a 100Mhz 486) would be regarded as archaic junk today. 

All of the examples are systems-level inverse problems, involving simultaneous equations and multiple dimensions, as this is the way problems are apprehended in science—not the simple explicit structure of spreadsheets and ordinary programming languages. Yet the complexity of programming, as illustrated above, is hardly more than what it takes to use a spreadsheet program.  Compare this to what it is like to program in mainstream languages like Perl, Java, or C++, and ask yourself how one may go about programming these applications in those languages!

The major lever for rapid prototyping is that in each holon of a program structure, the solver is interchangeable with other solvers in its class. Thus a user can experiment with solvers in the way he experiments with the multiple search-engines on the web.  After all, that is what implicit solvers are—search engines. This facility dramatically reduces the time it takes to prototype a new application.

Meta Calculus in Education

Although the level of problem solving is mathematically advanced, beyond the level learned in normal undergraduate engineering education, there is nothing in the use of this technology that requires that level of education. In point of fact, the only mathematical skill required is the ability to pose problems in algebra and trigonometry. A knowledge of calculus methodology is not required.

What this means is that meta calculus could be used as a tool for learning science without the necessity of understanding how problems are solved mathematically.  An approach stressing concepts of mathematical science, divorced from the methodology of math, could be undertaken at a very early age, in which meta calculus languages could be employed in ways that motivate understanding of science without requiring mathematical rigor that tends to demotivate students. This non-rigourous approach is normally taken in college physics courses for liberal arts students who have not been exposed to calculus or analytic geometry. It is also the approach taken in most high-school physics courses.

An example of how this might be used to motivate students to become enthusiastic about science would be to develop an online encyclopedia of pedagogical demo problems containing thousands of small application problems such as these, representing every branch of science and engineering. Each application would be modeled step by step from first principles, to illustrate and promote the understanding of science, and dynamic graphics would be used to illustrate the dynamic stages of the solution processes (like the pilot-ejection profile problem and the chemical kinetics problem), thereby illustrating the effect of hidden mathematical algorithms (like in online games), without immersing the student in the mathematical methodology. Mathematical methodology would thereby be relegated to electives for students interested primarily in mathematics or computer science. 

The key difference in the way math and science is taught today is that math methodology would pushed into the background and application science methodology brought to the foreground, rather than vice-versa which characterizes the pedagogy of the past. It is well known that online games are addictive like narcotics. This is certainly exploited by game producers. Why not make science education addictive in the same way by building exotic games from the pedagogy of science? In this way detailed understanding of science could be imparted to children in the primary grades and such understanding would become pervasive in the society as a whole within three generations.