Static/dynamic control for optimizing a useful objective

[0002] 1. Technical Field

[0003] The present invention relates to a method, computer system, and program product for optimizing a useful objective function with respect to a dynamic parameter vector and a static parameter vector, and more particularly, to optimizing a spacecraft objective such as a spacecraft trajectory objective and a spacecraft design objective.

[0004] 2. Related Art

[0005] A goal of spacecraft trajectory optimization is to obtain the best trajectory. Spacecraft trajectory optimization is a dynamic problem. A dynamic problem involves a sequence of decisions made over time. For example, the thrust direction for a spacecraft engine must be selected on day one of a space flight. Then a new thrust direction must be selected for day two, and so on. Each previous decision affects all future decisions. A complete set of decisions for the total flight time is necessary to define a trajectory.

[0006] A goal of spacecraft design optimization is to obtain the best design. Spacecraft design optimization is a not a dynamic problem. Design decisions are not linked to a time sequence. For example, the decisions of how large a solar array to build, how much fuel to carry, and when to launch a spacecraft are static or non-dynamic decisions.

[0007] The overall performance of a spacecraft is a function of both the spacecraft trajectory and the spacecraft design. Further, the spacecraft design will impact the spacecraft trajectory and vice versa. Since spacecraft design and spacecraft trajectory are not independent, the most efficient spacecraft design and the most efficient spacecraft trajectory cannot be determined independently.

[0008] Current methods for optimization are constructed to solve either dynamic problems or static problems. Related art optimization methods cannot solve problems which have both a static part and a dynamic part without approximation. Related art methods for optimization can be divided into two categories. The first category is non-dynamic or static optimization methods known as “parameter optimization”. The second category is dynamic optimization methods known as “optimal control.”

[0009] Parameter optimization was originally developed to solve static problems. Parameter optimization is directly applicable to static problems, such as spacecraft design. However, parameter optimization is not easily or accurately applicable to dynamic problems such as space flight trajectories. Parameter optimization can only be applied to dynamic problems when the true dynamic nature of the problem is approximated or removed entirely. Since parameter optimization cannot solve general dynamic problems without approximation, parameter optimization cannot solve the combined static and dynamic problem of spacecraft design and trajectory without approximation.

[0010] Optimal control was developed for exclusively dynamic problems. Optimal control is applicable to dynamic problems such as spacecraft trajectory optimization. However, optimal control is not suited to static problems such as spacecraft design.

[0011] The related art uses either parameter optimization or optimal control to design deep space missions. The related art cannot achieve the combined benefits of both parameter optimization and optimal control simultaneously. The following two subsections describe the parameter optimization and the optimal control methods presently in use.

[0012] The related art applies parameter optimization to the dynamic problem of spacecraft trajectory optimization by making several significant approximations. The physics of the trajectory problem is approximated by instantaneous impulses followed by ballistic coasts. The instantaneous impulses are intended to represent the effect of thrust on a spacecraft. A sequence of impulses and coasts can roughly approximate a continuously running engine. The ballistic coasts usually only account for the gravitational influence of a single central object (usually the sun.) All long range planetary gravitational effects are ignored. Since planetary gravitation is ignored, the ballistic coasts are perfect conic sections whose shape is always time independent. The time independence of the shape is a necessary assumption for the non-dynamic parameter optimization method.

[0013] Of particular interest is the ion propulsion engine. Ion propulsion engines are far more efficient than ordinary traditional chemical engines. Ion propulsion engines are expected to replace many applications of chemical engines in the near future. Ion propulsion engines have significantly different operating characteristics than chemical engines. In particular, ion engines typically operate continuously for days or even years at low thrust intensities. The impulse/ballistic coast approximation required by parameter optimization is a particularly poor approximation when applied to ion engines.

[0014] The parameter optimization method used in the related art requires prespecification of the sequence of planets that the spacecraft will pass by closely (“planetary flyby sequence”). Thus, the planetary flyby sequence is not susceptible to parameter optimization. Individual planetary flybys are prespecified by fixed constraints. Prespecification of flybys greatly reduces the likelihood of discovering the most efficient flyby sequence.

[0015] Parameter optimization typically requires the physics of planetary flybys to be approximated as a collision at a single point in space. The spacecraft trajectory is propagated to the center of the flyby body without accounting for the gravity of the body. Flybys are then modeled as an instantaneous change in velocity without a change in the spacecraft position. With parameter optimization, the true extended spatial and dynamic nature of flybys is not correctly represented. This approximation significantly reduces the precision of parameter optimization solutions.

[0016] Related art parameter optimization methods do not optimize spacecraft thrust sequences directly. Instead, related art methods optimize instantaneous changes to spacecraft velocity. Optimizing velocity changes neglects or only approximates the effect of the dynamic nature of the spacecraft mass. The spacecraft mass is a decreasing function of time because fuel is burned. Neglecting this fact reduces the precision of parameter optimization solutions.

[0017] The related art uses an optimal control method to calculate or fine tune navigation paths or trajectories for spacecraft. The optimal control method is known as the calculus of variations (“COV.”) The COV method is not capable of calculating optimal trajectories for spacecraft from scratch. An efficient trajectory and/or a prespecified flyby sequence must be supplied as input. The input trajectory, or trajectory associated with a prescribed flyby sequence, is not derived from a precise optimization procedure and is therefore not necessarily optimal. Similarly, the prespecified flyby sequence is not necessarily optimal.

[0018] Related art COV methods typically make several approximations similar to the approximations made by related art parameter optimization methods. For example, thrust is approximated as a sequence of impulses, spacecraft propagation only takes into account the Sun's gravity, and planetary flybys are treated as an instantaneous change in spacecraft velocity. As a result, the COV method is limited in precision in the same way that the parameter optimization method is limited in precision.

[0019] The COV method also suffers from extreme sensitivity to flyby parameters. The extreme sensitivity of COV results in a substantial reduction in the improvement that COV can potentially achieve when flybys are involved. The sensitivity problem limits the COV method to consider trajectories with only a small number of flybys.

[0020] The main advantage of the COV method is that it is a dynamic method. The dynamic aspect of the COV method could, in theory, permit the correct representation of the dynamics of space flight. Unfortunately, the COV method is not robust enough to solve the trajectory problem without relying on significant approximations or requiring the prespecification of the flyby sequence.

[0021] Application of the preexisting general optimization methods, known as neural networks, genetic algorithms, and simulated annealing to static/dynamic problems, may produce inferior results that do not converge to optimal solutions. These optimization methods are heuristic or approximate methods involving various types of random searches to find good solutions. Heuristic methods often fail to find the best or optimal solutions.

[0022] An accurate method is needed to simultaneously optimize both dynamic and static variables collaboratively, such as by simultaneously taking into account the effect of navigation on spacecraft design, and vice-versa.

[0023] The present invention of static/dynamic control (“SDC”) provides a method for optimizing a useful objective, comprising the steps of:

[0024] providing a computer system having:

[0025] a memory device;

[0026] a code located on the memory device;

[0027] a processor for executing the code; and

[0028] a useful objective function J coupled to the code, wherein the code includes a static/dynamic control program for optimizing J with respect to a dynamic control vector v(t) and a static parameter vector w, and wherein t denotes time;

[0029] providing input data for the static/dynamic control program, wherein the input data is coupled to the code;

[0030] optimizing J with respect to v(t) and w, by executing the static/dynamic control program and using the input data; and

[0031] outputting a computed result from the optimizing step to an output device within the computer system.

[0032] The present invention provides a first computer system for optimizing a useful objective, comprising:

[0033] a memory device,

[0034] a computer code located on the memory device;

[0035] a computer processor for executing the computer code;

[0036] a useful objective function J coupled to the computer code, wherein the computer code includes a static/dynamic control program for performing an optimization of J with respect to a dynamic control vector v(t) and a static parameter vector w, and wherein t denotes time; input data for the static/dynamic control program; and an output device for receiving a result from the static/dynamic control program.

[0037] The present invention provides a second computer system for optimizing a useful objective, comprising:

[0038] a memory means for storing data on a memory device;

[0039] a computer code, located on the memory device, for optimizing a useful objective function J with respect to a dynamic control vector v(t) and a static parameter vector w, wherein t denotes time;

[0040] a processor means for executing the computer code;

[0041] an input means for transferring input data to the static/dynamic control program; and

[0042] an output means for receiving a result from the static/dynamic control program.

[0043] The present invention provides a computer program product, comprising:

[0044] a recordable medium; and

[0045] a computer code recorded on the recordable medium. wherein the computer code includes a static/dynamic control program.

[0046] By optimizing a useful objective, the present invention produces a useful, concrete, and tangible result. In particular, the useful objective may include a spacecraft objective, a groundwater decontamination objective, or a chemical reactor objective.

[0047] The present invention has the advantage of solving problems having both a static and a dynamic aspect. In the related art, optimization algorithms are designed for either static problems or dynamic problems, but not both.

[0048] The present invention has the advantage of solving a generalized state equation. The state equation include such equations as Newton's equations of motion, general relativistic equations of motion, Darcy's law of potential flow for groundwater movement, computational fluid dynamics equations for air flow over an aircraft, and transient equations for the transport of neutrons, radiation, heat, electric current, etc.

[0049] The present invention has the advantage of accommodating generalized physical constraints such as requiring: a minimum distance between a spacecraft and the sun so as to protect the spacecraft from overheating, a maximum feasible thrust magnitude for a spacecraft, or a maximum allowed temperature during operation of a chemical reactor in a chemical processing system.

[0050] The present invention has the advantage of allowing the initial or the terminal condition of the state of the system to be a variable. Related art optimal control methods require the initial state of the system to be given and fixed.

[0051] The present invention has the advantage of including an algorithm having two formulations: a period formulation and a fully continuous formulation. With the period formulation, the dynamic control is permitted to be a discontinuous, parametric function of time. The usefulness of the period formulation is illustrated by the capability of changing engine throttle once every fixed number of days, such as once per ten days. With the fully continuous formulation, the dynamic control is permitted to be non-parametric and continuous in time. Illustrative of the usefulness of the fully continuous formulation is the capability of satisfying a generally recognized need to optimize solar sail and light propulsion spacecraft, which are expected to be deployed in the next few years.

[0052] The present invention has the advantage of solving spacecraft optimization problems without approximation to the physics or dynamics involved. Among the effects susceptible of inclusion by SDC are general relativistic corrections, radiation pressure, continuous thrust, planetary harmonics, dynamic interaction between spacecraft and planetary flybys, gravitational influence of any number of solar system objects, and planetary drag on spacecraft traversing a planetary atmosphere. Related art methods require simplifying assumptions for compatibility with related art algorithms.

[0053] The present invention has the advantage of satisfying a generally recognized and existing need to optimize both ion propulsion spacecraft and combination chemical and ion propulsion spacecraft designs and trajectories. Related art methods perform particularly poorly for ion propulsion and general low thrust spacecraft.

[0054] The present invention has the advantage of enabling mathematical derivatives to be computed exactly. Exact derivatives result in a more robust process and faster convergence. The related art typically uses approximate derivatives.

[0055] The present invention has the advantage of employing a fully second order optimization algorithm which cannot converge to false minima or saddle points. Both the COV and parameter optimization methods currently in use can converge to false minima or saddle points.

[0056] The present invention has the advantage of converging only to optimal solutions which satisfy both the necessary and sufficient conditions of optimality. Existing heuristic optimization methods (neural networks, genetic algorithms, simulated annealing) often fail to converge to optimal solutions, since they are not designed to find optimal solutions.

[0057] FIG. 1 depicts a flow chart of the SDC method, in accordance with a preferred embodiment of the present invention.

[0058] FIG. 2 depicts the SDC algorithm of FIG. 1 for the period formulation.

[0059] FIG. 3 depicts an example of a single component of the dynamic control vector v(t) that is constant within each of several time periods, in accordance with the period formulation of FIG. 2.

[0060] FIG. 4 depicts the period formulation of FIG. 2 applied to an example of spacecraft design and trajectory optimization.

[0061] FIG. 5 depicts the first-iteration spacecraft trajectory for the example in FIG. 4 with the initial iteration values of: spacecraft travel for 3.25 years beginning Jan. 1, 2014 with zero spacecraft thrust for the entire flight.

[0062] FIG. 6 depicts the optimal converged spacecraft trajectory for the example of FIGS. 4-5, using nuclear electric propulsion (NEP).

[0063] FIG. 7 depicts the time dependence of the converged optimal thrust sequence and velocity magnitude for the example of FIG. 6.

[0064] FIG. 8 depicts an assumed maximum thrust verses distance to the sun for the example of FIG. 9, using solar electric propulsion (SEP).

[0065] FIG. 9 depicts the optimal converged spacecraft trajectory for the example of FIGS. 4-5, using solar electric propulsion (SEP).

[0066] FIG. 10 depicts the time dependence of the converged optimal thrust sequence and velocity magnitude for the example of FIG. 9.

[0067] FIG. 11 depicts the SDC algorithm of FIG. 1 for the fully continuous formulation.

[0068] FIG. 12 depicts a flow chart of a computer system showing the data coupling associated with computer implementation of the SDC method, in accordance with a preferred embodiment of the present invention.

[0069] This section presents the form and theoretical basis for a novel process called the Static/Dynamic Control (SDC) process which solves useful physical problems. The general physical problem's form arid the SDC process which solves these problems was initially invented to address spacecraft design and trajectory optimization problems. However, the SDC process is applicable to a broad class of physical problems in addition to spacecraft design and trajectory optimization, including groundwater decontamination system design and operation, oil extraction system design and operation, chemical reactor design and operation, economic decision modeling, aircraft design and operation, automobile design, reservoir system design and operation, agricultural crop yield maximization, agricultural pest minimization, and manufacturing process design.

[0070] The SDC process optimizes a useful objective. A useful objective is a state or condition expressing a magnitude or quantity associated with a practical application. An example of a useful objective is the volumetric size of an automobile engine. Another example of a useful objective is the mass of fuel utilized in a spacecraft mission. A useful objective is mathematically represented by a useful objective function J, such as the function J defined in Equation (2) infra. J is a function of optimization parameters. Optimizing the useful objective, or optimizing the useful objective function J, is accomplished by either maximizing or minimizing J with respect to the optimization parameters. A value of J is defined as the computed value of J when J is evaluated at particular values of the optimization parameters.

[0071] The SDC process includes an algorithm having two formulations: a period formulation and a fully continuous formulation. With the period formulation, the dynamic optimization variables are permitted to be discontinuous in time, and are alternatively permitted to be continuous in time. With the fully continuous formulation, the functions serving as the dynamic optimization variables are represented only as continuous in time.

[0072] Utilization of the SDC process includes the following steps. First, a physical problem is posed in a form which is analogous to the SDC period formulations of subsections 4 and 6, or the fully continuous formulation of Appendix B, subsection 2. Second, the problem is solved using the SDC algorithm steps specified in subsection 5 for the period formulation or in Appendix B, subsection 3 for the fully continuous formulation. Third, the physical design produced by the SDC algorithm is implemented.

[0073] Part of the SDC process is the SDC algorithm. The SDC algorithm addresses situations when it is desirable to optimize a set of time-varying variables, and a set of fixed-in-time variables simultaneously with the added complexity of letting the state at either the start or the end of time be a function of control variables. In addition, the algorithm described here is capable of optimizing the time-duration and start-time for general physical problems. When the SDC process is applied to spacecraft, the spacecraft design, launch date, launch energy, flight-time and trajectory can be optimized simultaneously.

[0074] While the SDC method can be implemented analytically for simple problems, practical problems typically require computer implementation. The data coupling associated with computer implementation for computer system 600 is illustrated in FIG. 12. Input data 620 is coupled to SDC computer code 610. An optimization of a useful objective function J, which is defined and discussed in subsection 3 infra, is accomplished by execution of the SDC computer code 610, using the algorithm of the SDC method. The SDC computer code 610 may be located within a first portion 640 of a memory device 630 of the computer system 600. The computer system 600 includes a processor 670 for executing the SDC computer code 610. As a result of the execution of the SDC computer code 610, output data 680 is generated. The memory device 630 may include, inter alia, read-only memory (ROM), random-access memory (RAM), or both. The computer code 610 could be hard wired into ROM or loaded into RAM, of the memory device 630, as needed. The coupling of the input data 620 to the SDC computer code 610 can be accomplished many different ways. One way is to have some or all of the input data 620 hard wired into the SDC computer code 610. Another way is to have the input data 620 stored a second portion 650 of the memory device 630. The coupling of the input data 620 to the SDC computer code 610 can also be accomplished by transfer of the input data 620 from an input device 615 to the SDC computer code 610. The input device 615 could be any input device, such as a computer keyboard, a removable or nonremoveable hard disk drive, a floppy disk drive, a Bernoulli drive, or a magnetic tape drive. The output data 680 may include any data available to the SDC computer code 610, such as a computed result by the SDC computer code 610. The output data 680 could be transferred from the computer code 610 to any output device 690, such as a printer, a plotter, a removable or nonremovable hard disk drive, a floppy disk drive, a Bernoulli drive, a magnetic tape drive, or the screen of a computer monitor. The SDC computer code 610 could be recorded on any recordable medium, such as the memory device 630, a disk, or a magnetic tape.

[0075] The various components, such as the memory device 630, processor 670, input device 615, and the output device 690, within the computer system 600 may each be physically located at any distance with respect to one another. For example, the memory device 630 and processor 670 may be located in Chicago, the input device 615 may be located in Denver, and output devices 690 may be located in both Denver and Boston.

[0076] The useful objective function J may be coupled to the SDC computer code 610 in any practical manner. For example, the useful objective function J may be hard wired into the SDC computer code 610. As another example, a plurality of useful objective functions J may be included in SDC computer code 610, with the user of the SDC computer code 610 selecting his or her desired J by means of user input data 620. A third example is having the useful objective function J formed by a code module that is external to the SDC computer code 610 and is executable by the SDC computer code 610. A fourth example is having the SDC computer code 610 choose one of several available useful objective functions J as a function of the result of an auxiliary calculation by the SDC computer code 610. For example with a spacecraft application involving Mars as a flyby, the SDC computer code 610 could choose to minimize the amount of fuel spent, or to maximize the time during which the spacecraft is within a given distance from the surface of Mars, as a function of the computed distance between Mars and the Sun when the spacecraft is within the gravitational environment of Mars. Numerical parameters associated with a useful objective function J may be coupled to the computer system 600 in various ways, such as by being hard wired into the SDC computer code 610, stored on an auxiliary storage device such as a hard disk, or supplied by user input such as by user keyboard input.

[0077] An applicable computer system 600 includes any mainframe computer or desktop computer, inasmuch as the computations for many applications can be performed in a practical amount of time with a personal computer. For the Example 2 discussed infra in subsection 3, part C, the SDC calculations comprising 1941 iterations were performed in only 128.5 minutes on a SUN Ultra Sparc workstation.

[0078] Subsection 2 presents the notation conventions used in the description and theoretical basis of the SDC algorithm. A formulation of SDC applied to general engineering problems is given in subsection 3, part A. Subsection 3, part B presents an example of how the period formulation of SDC is applied to spacecraft design and trajectory problems. Two numerical examples of spacecraft optimization are given in subsection 3, part C. A detailed description of the SDC algorithm for the period formulation, applicable to general engineering problems. is presented in subsections 4 through 6. The theoretical basis of the SDC period formulation is presented in subsections 7 and 8. Appendix B describes the fully continuous formulation of the SDC algorithm. The algorithm steps and theoretical basis in Appendix B apply to general engineering problems which have time-continuous controls (i.e. period functions or dynamic limitations are unnecessary.)

[0079] The following notation convention is used. A general variable or function may be either upper or lower case letter. A variable or function may be scaler, vector, matrix, or rank 3 valued. The dimensions of a variable or function are determined by the definition and context. A complete list of definitions is included as appendix A.

[0080] The following general rules apply when a variable or function is referenced. Round brackets always enclose arguments for functions. For example, g(2) is the value of the function g evaluated with the argument 2. Square brackets enclose element indices for variables or functions of rank greater than zero. For example T[j] is the j^thcomponent of the vector T. Square and round brackets can be combined. For example, T(x, v, w, t)[i] is the i^thcomponent of the vector fiction T evaluated using the arguments x, v, w, and t.

[0081] Superscripts denote ordinary powers for rank zero objects. A superscript “t” on rank one and two objects indicates the transpose operation. Row and column vectors are distinguished. A vector referenced without a superscript t is assumed to be a column vector. A superscript of “−1” on rank two objects indicates the matrix inverse operation.

[0082] Subscripts involving x, v, u, w and/or t imply the object is a function and has been differentiated. For example $a_{\mathrm{ux}}=\frac{{\partial }^{2}a}{\partial u\partial x}.$

[0083] The order of a differentiated function's elemental indices (in square brackets) is first undifferentiated function indices followed by indices resulting from differentiation. For example suppose T is a vector function. A rank three object is obtained by differentiating T with respect to the two vectors w and x. A single element of the rank three object is: $\begin{array}{cc}{T}_{\mathrm{wx}}\left(x,v,w,t\right)\left[i,j,k\right]=\frac{{\partial }^{2}T\left[i\right]}{\partial w\left[j\right]\partial x\left[k\right]}.& \left(1\right)\end{array}$

[0084] An element index of “:” indicates this index is free. For example, a rank two object (matrix) based on Equation (1) is

[0085] T_wx(x, v, w , t)[i, :, :]=∇_w(∇_xT[i]))

[0086] A specific part of an object can be specified using the colon notation. For example, x(t)[4:6] is a vector with three components consisting of the fourth, fifth, and sixth components of the vector x(t). Similarly, T_wx(x, v, w, t)[i, 1:2, 1:3] is a two by three matrix consisting of the upper two by three block of the matrix T_wx(x, v, w, t)[i, :, :].

[0087] Ordinary scaler, vector, and matrix multiplication is implied when rank zero, one and two objects are multiplied. For example, suppose a is a scaler function and T is a vector function of vectors x, u, and w, then the [i, j] component of the product a_uxT_wis interpreted as: $\mathrm{product}\left[i,j\right]=\sum _{k=1}^n\frac{{\partial }^{2}a}{\partial u\left[i\right]\partial x\left[k\right]}·\frac{\partial T\left[k\right]}{\partial w\left[j\right]}$

[0088] Similarly, the [i, j] component of the product a_uxT_w^tis interpreted as: $\mathrm{product}\left[i,j\right]=\sum _{k=1}^n\frac{{\partial }^{2}a}{\partial u\left[i\right]\partial x\left[k\right]}·\frac{\partial T\left[j\right]}{\partial w\left[k\right]}$

[0089] Multiplication involving rank three objects is indicated with explicit summations to avoid ambigtuity.

[0090] The symbol R represents the set of all real numbers. The symbol R is used to indicate the dimension of the value of all functions and variables. For example, the statement “uεR^m” indicates u is a vector with m components or is m-dimensional. The statement “βεR^m×n” indicates β is a matrix with m rows and n columns. Expressions of the form 0^m×nrepresent a matrix of all zeros with m rows and n columns.

[0091] Integer subscripts i, j, or 1, 2, etc. arc interpreted as part of a variable's name. Integer subscripts do not indicate differentiation. Subscripts of this type only appear as part of the names of the variables u_i, t_i, α_i, β_i, and γ_i(to be defined later.) For instance, β₂is a matrix variable εR^m×n. The j^thcolumn of β₂is represented by β₂[:, j] which is εR^m. Subscripts of integer type are not used as part of function names.

[0092] The symbol “≐” is read “is defined to be” to distinguish between definitions and ordinary “=” assignment or equality statements. A “bar” over a function (for example {overscore (J)}_uu) indicates the function is evaluated at a given or “nominal” point. A function with a bar over it is at most only a function of time because the dependence on all other variables has been eliminated. The “Hessian” of a scaler function is defined to be the second derivative of the scaler function with respect to a vector. In general, a Hessian is a rank two object. The expression “H.O.T.” stands for “Higher Order Terms.” A complete list of variable and function definitions is included as Appendix A.

[0093] Either the period formulation or the fully continuous formulation of the SDC method of the present invention can be applied to any physical problem in which static and dynamic parameters are to be simultaneously optimized. Examples of engineering problems to which SDC can be applied are groundwater decontamination system design and operation, oil exploration system design and operation, spacecraft design and operation, chemical reactor design and operation, economic decision modeling, and aircraft design and operation. In particular, the SDC process can be applied to any engineering problem which can be described by three general equations, such as Equations (2), (3), and (4) infra.

[0094] FIG. 1 shows a flow diagram of the SDC method for both the period formulation and the fully continuous formulation. FIG. 1 illustrates how SDC is applied to a general engineering problem such as a static design and dynamic control problem. A useful objective 400 for the period formulation is represented by the following useful objective function J: $\begin{array}{cc}J={\int }_{t_o}^{t_N}F\left(x\left(t\right),v\left(t\right),w,t\right)dt+\sum _{i=1}^{N}G\left(x\left(t_i\right),v\left(t_i\right),w,t_i,i\right),& \left(2\right)\end{array}$

[0095] or a mathematical equivalent. Both t₀<t_Nand t₀>t_Nare within the scope of the SDC period formulation. Note that the fully continuous formulation permits an analogous scope of time ordering (t₀<t_fand t₀>t_f, as discussed in Appendix B). Any mathematical equivalent can be used in place of Equation (2). The analogous expression for J in the fully continuous formulation is Equation (182) in Appendix B. The goal of SDC is to optimize J by choosing the optimal or “best” dynamic control vector v(t) and optimal or “best” static parameter vector w. The useful objective J can be either minimized or maximized in value. The functions F and G in Equation (2) (or Equation (182) in the fully continuous formulation) are selected by the user of the invention to best represent the design and control objectives, Examples of optimizing useful objectives include minimizing spacecraft fuel usage, minimizing the contamination remaining in the ground for environmental decontamination problems, or maximizing useful products from a chemical reactor.

[0096] In the period formulation of SDC, the functions F and G can depend on the state vector x(t), the dynamic control vector v(t), the static parameter vector w, and time t. In the fully continuous formulation of SDC, F has the same dependencies as in the period formulation. and C can depend on w, x(t_f), v(t_f), and t_f. The time t_fis a “final,” or terminal, value of t if t₀<t_f. However, t_fis an “initial” value of t if t₀>t_f.

[0097] The state vector x(t) represents the physical state of the engineering system at time t. The state vector encompasses variables not under the direct control of the engineer. Examples of x(t) include temperature, velocity, and chemical disequilibrium. In general, the engineer has indirect control over state variables x(t), such as temperature. For example, the engineer can change temperature by adjusting a thermostat. However, temperature can not be changed instantaneotusly, since temperature varies continuously according to the equations of transient heat transfer. The dynamic control vector v(t) represents the dynamic parameters which are under direct control by the engineer. Examples of v(t) include a thermostat setting, rocket thrust magnitude, and chemical inflow rate in a chemical reactor. The dynamic control vector v(t) represents decisions which may change continuously or discontinuously in time. Like v(t), the static parameter vector w represents parameters which are under direct control of the engineer. Unlike v(t), however, the static parameter vector w is associated with engineering design decisions which are not functions of time. Examples of uw include the size of a tank to build for a chemical reactor, the size of a pipe to use, or the launch date for a spacecraft. In general, the user of the invention defines v(t) and w to be the dynamic and static parameters for which optimal values are desired.

[0098] The second equation of the SDC process, in either the period formulation or fully continuous formulation, is the system state equation 401. The state equation includes the general form, $\begin{array}{cc}\frac{dx\left(t\right)}{dt}=T\left(x\left(t\right),v\left(t\right),w,t\right)x\left(t=t_0\right)=\Gamma \left(w\right),& \left(3\right)\end{array}$

[0099] or a mathematical equivalent. The vector function T is selected by the user of the invention to best represent the time evolution of the state vector x(t) under the influence of the current state, the current dynamic control vector v(t), and the static parameter vector w. In both the period formulation and the fully continuous formulation, the vector function T may be a discontinuous function of time. The vector function Γ(w) represents an initial condition on x(t) if t₀<t_fand a terminal condition on x(t) if t₀>t_f. The condition of the system can be optimized because it is allowed to be a function of the static parameter vector w. Examples of system state equations include: the heat equation, chemical equilibrium equations, groundwater contamination transport equations, dynamic equilibrium equations, equations describing the flow of electricity, and Newton's laws of motion under the influence of forces.

[0100] The third equation of the SDC process is the constraint equation 402. Use of the constraint equation 402 is optional, since the SDC formulation works with or without the constraints 402. The constraint equation 402 represents any physical or practical engineering constraints that may exist for a particular problem. Such engineering constraints are written in the general form,

[0101] L(x(t), v(t), w, t)≧0 and/or K(x(t), v(t), w, t)=0,  (4)

[0102] or a mathematical equivalent. Equation (4) applies to either the period formulation or the fully continuous formulation. The linear or nonlinear vector functions L and K are selected by the user of the invention to represent practical or physical constraints on the engineering problem. Examples of constraints include a maximum temperature for a chemical reactor, a maximum allowed level of residual contamination for environmental decontamination problems, or a minimum allowed distance between the sun and a spacecraft to avoid spacecraft overheating.

[0103] The fourth equation is used only by the period formulation of the SDC process and is the dynamic limitations equation 403. The dynamic limitations equation 403 is optional, since the SDC formulation works with or without dynamic limitations. The fully continuous formulation of the SDC process does not use the dynamic limitations equation 403. The dynamic limitation equation 403 represents any physical or practical engineering constraints on the time evolution of the dynamic control vector v(t) that may exist for a particular problem. The period formulation of the SDC process requires that the dynamic limitations be written in the general form, $\begin{array}{cc}v\left(t\right)=\left|\begin{array}{cc}f\left(u_1,w,t,1\right)& \mathrm{for}t=t_0\mathrm{to}{t}_1\\ f\left(u_2,w,t,2\right)& \mathrm{for}t=t_1\mathrm{to}{t}_2\\ ⋮& ⋮\\ f\left(u_N,w,t,N\right)& \mathrm{for}t=t_{N-1}\mathrm{to}{t}_N.\end{array}\right|& \left(5\right)\end{array}$

[0104] The vector functions 𝑓(u_i, w, t, i) are selected by the user of the invention to properly represent any mechanical or practical limitations on the time evolution of v(t). The number of periods N may be chosen arbitrarily. The functions 𝑓 are parameterized by the dynamic parameter vector u_i, the static parameter vector w, and time t. The functions 𝑓 can be used to effectively limit or constrain the SDC algorithm to consider only solutions v(t) which are of the form of Equation (5). The time intervals t_ito t_i+1are called “periods.” The dynamic parameter vector u_iis constant within each period i, i=1, 2, . . . , N. For example, the simplest useful set of functions 𝑓 is 𝑓(u_i, w, t, i)=u_i. The dynamic control vector v(t) may be optimized such that v(t) is constant over each period, allowing changes only at period interfaces t_i. Alternatively, v(t) may be subject to a dynamic limitation that allows v(t) to vary within each period, either continuously or discontinuously. An example of a single component of v(t) for the special case when v(t) is required to be constant within each period is plotted in FIG. 3. The example in FIG. 3 corresponds to a situation where it is impractical or infeasible to adjust a dynamic control continuously in time. An example is when v(t) is spacecraft thrust and commands to change spacecraft thrust can only be transmitted once every day. An example of incorporating constraints 402 into the useful objective 400 and incorporating dynamic limitations 403 into “period” functions 𝑓 is given in the subsection 3, part C.

[0105] Once the user of the invention has selected a useful objective function 400, a state equation 401, (optionally) a constraint equation 402, and a dynamic limitation equation 403, the static/dynamic algorithm 404 can be used. The static/dynamic (SDC) algorithm 404 generates the optimal values for v(t) and w such that the useful objective 400 is optimized while obeying the state equation 401, and any constraints 402 and dynamic limitations 403. The optimum values of v(t) and w are used in step 405 for obtaining an optimum static engineering design and/or an optimum dynamic control.

[0106] FIG. 2 is a flow chart describing the period formulation of the SDC algorithm 404 in greater detail. The SDC algorithm for the period formulation 404 is an iterative procedure which includes an inner loop terminating in an adequacy of improvement test 411, and an outer loop terminating in a convergence test 412. The step 406 in the period formulation of the SDC algorithm 404 incorporates the constraints 402 into the useful objective 400 and/or incorporates dynamic limitations 403 into “period” functions 𝑓 (see FIG. 1 for useful objective 400, constraints 402, dynamic limitations 403, and SDC algorithm 404). Note that incorporation of dynamic limitations 403 into period functions 𝑓 is implemented only if the SDC application includes a dynamic limitation.

[0107] The step 407 in FIG. 2 of the period formulation of the SDC algorithm 404 provides arbitrary initial approximations for v(t) and w, and then uses the initial approximations of v(t) and w and the state equation 401 to compute x(t). The state x(t) is obtained by integrating Equation (3) using any known method for integrating a differential equation subject to a single time condition. Once x(t) is obtained, the value of the useful objective function J is calculated using the useful objective 400 as represented by Equation (2).

[0108] The step 408 of the period formulation of the SDC algorithm 404 initializes program parameters such as P, Q, R, S, Y, W, and c as discussed in subsection 5 infra.

[0109] The step 409 in the SDC algorithm 404 integrates the static/dynamic control equations, i.e., Equations (32)-(40) infra, to obtain an update law. An update law is a prescription for updating values of v(t), w, and x(t), such as equations of update (50)-(52), which are discussed in subsection 5 infra. Obtaining an update law includes calculating parameters, such as α_i, β_i, and γ_i(i=1, 2, . . . , N), for use in equations of update.

[0110] The step 410 in FIG. 2 of the period formulation of the SDC algorithm 404 uses the update law obtained from the step 409 and the state equation 401 to compute updated values for v(t), w, and x(t). Once updated values for v(t), w, and x(t) are obtained, an updated value of J is calculated using the useful objective 400 as represented by Equation (2).

[0111] The step 411 in FIG. 2 of the SDC algorithm 404 tests for adequate improvement in the values of J. An improvement is adequate if the updated value of J has changed sufficiently from the value of J at the previous iteration. The adequacy of improvement test is mathematically articulated in Equation (53) infra. If the improvement is adequate, then the step 412, the convergence test as stated in subsection 5 infra, is executed. If the improvement is not adequate, then step 410 is repeated with an increasingly damped update law. The update laws are damped by decreasing the value of the parameter ε as described in subsection 5 infra.

[0112] The step 412 in FIG. 2 of the SDC algorithm 404 tests for convergence to the optimal values of v(t) and w. If convergence is obtained, then the procedure stops. If convergence is not obtained, then the values of v(t), w, and x(t) are replaced with the improved values obtained in the steps 410 and 411, followed by a return to the step 408.

[0113] By optimizing the useful objective function J, the SDC method of the present invention produces a useful, concrete, and tangible result. In particular, the useful objective may include a spacecraft objective, such as a spacecraft trajectory objective and a spacecraft design objective. Optimizing a spacecraft trajectory objective provides spacecraft design and control of a spacecraft trajectory, including a capability of optimizing such parameters as flyby sequence, launch date, flight duration, and spacecraft, velocity. Optimizing a spacecraft design objective facilitates design of such parameters as engine efficiency, engine size, propellant tank size, and electric battery size.

[0114] A useful objective may also include an objective related to groundwater or soil decontamination. An example of a useful objective function J is a measure of the contamination remaining in the ground after a certain period of time. The SDC process can be used to design and control a decontamination system to minimize the contamination present in the ground. The static parameter vector w may include the location and size of groundwater extraction wells and clean water injection wells, and treatment plant size. The dynamic control vector v(t) may include dynamic pumping schedules for all wells. Other remediation system design parameters which can be optimized include the size and location of vapor extraction equipment, bioremediation equipment, horizontal wells, and slurry walls.

[0115] Another application involving a useful objective includes stabilizing vibrations in structures. For example, a useful objective function J includes a measure of the amplitude of the vibrations. The SDC process can be used to design and control structures to minimize the amplitude of destructive vibrations. The dynamic control vector v(t) may include the effect of mechanical actuators to damp vibrations. The static parameter vector w may include the size and location of various fixed elements of the structure such as beams and trusses.

[0116] Another useful objective includes maximizing the value of a portfolio. For example, a useful objective function J includes the dollar value of a portfolio after a specified period of time. The SDC process can be used to design trading laws to maximize return to risk ratios. The dynamic control vector v(t) may include the effect of decisions to buy or sell equities. The static parameter vector w may include the fraction of the portfolio to be exposed to different levels of risk.

[0117] Another application relating to a useful objective includes electric circuit design. For example, a useful objective function J includes the circuit response time. The SDC process can be used to design and operate an electric circuit to minimize circuit response tine. The dynamic control vector v(t) may include the parameters associated with the circuit controller. The static parameter vector w may include maximum voltage and current tolerances, and bandwidth for various components of the circuit.

[0118] Another application relating to a useful objective includes the design and operation of a chemical reactor. For example, a useful objective function J includes the weight of the useful chemical products. The SDC process can be used to design and operate a chemical reactor in order to maximize the weight of the useful chemical products produced in the reactor. The dynamic control vector v(t) may include the inflow rate of raw chemicals into the reactor, thermostat settings, and pressure levels. The static parameter vector w may include the size of the reactor's chambers and the length of time the reactor is operated.

[0119] Another application relating to a useful objective includes the design and operation of a water reservoir system. For example, a useful objective function J includes the amount of hydroelectric power which can be generated by the release of reservoir water. The SDC process can be used to design and operate a system of reservoirs in order to maximize the amount of hydroelectric power produced by the reservoir system. The dynamic control vector v(t) may include the flow rate of water into and out of the reservoirs in the system. The static parameter vector w may include the size of hydroelectric generators to build and the capacity or dam height for each reservoir in the system.

[0120] To enable the present invention to be better understood, this subsection describes an application of the present invention directed to optimizing spacecraft design and trajectory. The spacecraft design and trajectory problem is merely illustrative of the types of spacecraft and general engineering problems to which Static/Dynamic Control can be applied.

[0121] A manner of applying the period formulation of the SDC process to spacecraft design and trajectory optimization involves choosing a useful objective 400 represented as: $\begin{array}{cc}J={\int }_{w\left[1\right]}^{w\left[1\right]+w\left[2\right]}F\left(x\left(t\right),v\left(t\right),w,t\right)dt+\sum _{i=1}^{N}G\left(x\left(t_i\right),u_i,w,t_i,i\right),& \left(6\right)\end{array}$

[0122] where the scaler functions F and G are selected by the user of the invention to represent the spacecraft trajectory and design objectives. The system state vector x(t) contains at least six elements corresponding to the current spacecraft position and velocity. More state variables can be added if desired. An example of an additional state variable which may be appropriate is the current mass x(t)[7]=M_s(t) of the spacecraft. The dynamic control vector v(t) contains three components, one for each of the current x, y, and z components of thrust intensity or the extent and angular orientation of a solar sail. Defining v(t) in this way allows the spacecraft propulsion and trajectory to be optimized. Additional components of the dynamic control vector v(t) can be assigned to solar panel angle, spacecraft power consumption, etc. The static parameter vector w typically contains two components. The first, w[1], is the launch time or time t₀when the trajectory to be optimized is started. The second component, w[2], is the trajectory time duration or flight time. These definitions of w[1] and w[2] result in optimizing both launch date and flight time. Additional components of the static parameter vector w may include, engine performance parameters, solar panel size, launch energy, etc.

[0123] FIG. 4 shows the flow chart of FIG. 1 for an example of how SDC can be applied to spacecraft design and trajectory optimization. In the example illustrated by FIG. 4, the useful objective 413 is the mass of fuel used for the mission. Minimizing the mass of fuel used maximizes the payload of the spacecraft. The useful objective function J for this example can be written $\begin{array}{cc}J={\int }_{w\left[1\right]}^{w\left[1\right]+w\left[2\right]}\frac{d{M}_s\left(t\right)}{dt}dt.& \left(7\right)\end{array}$

[0124] In this case, $F=\frac{d{M}_s\left(t\right)}{dt}$

[0125] is the mass of fuel consumed per time and G is set to zero (i.e. G is not used.) The integral of F is the total mass of fuel consumed. Thus, by minimizing J, the total mass of fuel consumed is minimized.

[0126] The user of the invention may elect to limit the candidate spacecraft designs and trajectories by listing linear or nonlinear constraints 415 of the form of Equation (4). General examples of constraints 415 include a limit on how close a spacecraft is allowed to get to a planet or the sun, a maximum possible thrust intensity, a final destination, an intermediate destination, a maximum allowed flight time, an acceptable launch date range, a solar electric propulsion constraint, and a launch vehicle performance curve. For the example indicated by FIG. 4, the constraints 415 can be written,

[0127] x(w[1]+w[2])[1:3]−r_Saturn(w[1]+w[2])=0  (8)

[0128] w[1]−2001≧0  (9)

[0129] where r_Saturn(t) is the position of a target, such as Saturn, at time t and x(t)[1:3] (the first three components of x(t)) is the position of the spacecraft at time t. In this case, the function K=x(w[1]+w[2])[1:3]−r_Saturn(w[1]+w[2]) and the function L=w[1]−2001 in Equation (4).

[0130] In addition to providing a useful objective function 413 and constraints 415, a state equation 414 must be provided by the user of the invention in the form of Equation (3). The user of the invention must select appropriate vector functions T and Γ to properly represent the physics of the spacecraft problem. For the example indicated by FIG. 4, the user of the invention could define the state vector x(t) and the dynamic control vector v(t) as follows: $\begin{array}{cc}x=\left[\begin{array}{c}x\left[1\right]\\ x\left[2\right]\\ x\left[3\right]\\ x\left[4\right]\\ x\left[5\right]\\ x\left[6\right]\\ x\left[7\right]\end{array}\right]=\left[\begin{array}{c}x\mathrm{coordinate}\mathrm{of}\mathrm{spacecraft}\\ y\mathrm{coordinate}\mathrm{of}\mathrm{spacecraft}\\ z\mathrm{coordinate}\mathrm{of}\mathrm{spacecraft}\\ x\mathrm{velocity}\mathrm{of}\mathrm{spacecraft}\\ y\mathrm{velocity}\mathrm{of}\mathrm{spacecraft}\\ z\mathrm{velocity}\mathrm{of}\mathrm{spacecraft}\\ M_s\mathrm{mass}\mathrm{of}\mathrm{the}\mathrm{spacecraft}\end{array}\right]& \left(10\right)\\ v=\left[\begin{array}{c}v\left[1\right]\\ v\left[2\right]\\ v\left[3\right]\end{array}\right]=\left[\begin{array}{c}x\mathrm{component}\mathrm{of}\mathrm{thrust}\\ y\mathrm{component}\mathrm{of}\mathrm{thrust}\\ z\mathrm{component}\mathrm{of}\mathrm{thrust}\end{array}\right]& \left(11\right)\end{array}$

[0131] The state equation 414 for the example indicated by FIG. 4 is Newton's law of gravitation for m_pgravitating bodies. Thus the state equation 414 in the form of Equation (3) is $\begin{array}{cc}\frac{dx}{dt}=T\left(x,v,w,t\right)=\left[\begin{array}{c}x\left[4\right]\\ x\left[5\right]\\ x\left[6\right]\\ \frac{v\left[1\right]}{x\left[7\right]}+U_G\sum _{i=1}^{m_p}\frac{M\left(i\right)r\left(i\right)\left[1\right]}{{\Vert r\left(i\right)\Vert }^3}\\ \frac{v\left[2\right]}{x\left[7\right]}+U_G\sum _{i=1}^{m_p}\frac{M\left(i\right)r\left(i\right)\left[2\right]}{{\Vert r\left(i\right)\Vert }^3}\\ \frac{v\left[3\right]}{x\left[7\right]}+U_G\sum _{i=1}^{m_p}\frac{M\left(i\right)r\left(i\right)\left[3\right]}{{\Vert r\left(i\right)\Vert }^3}\\ \frac{\Vert v\Vert }{g{I}_{\mathrm{sp}}}\end{array}\right]=\left[\begin{array}{c}x\mathrm{velocity}\mathrm{of}\mathrm{spacecraft}\\ y\mathrm{velocity}\mathrm{of}\mathrm{spacecraft}\\ z\mathrm{velocity}\mathrm{of}\mathrm{spacecraft}\\ x\mathrm{acceleration}\mathrm{of}\mathrm{spacecraft}\\ y\mathrm{acceleration}\mathrm{of}\mathrm{spacecraft}\\ z\mathrm{acceleration}\mathrm{of}\mathrm{spacecraft}\\ \mathrm{rate}\mathrm{fuel}\mathrm{is}\mathrm{used}\end{array}\right],& \left(12\right)\end{array}$

[0132] where U_Gis the universal gravitational constant, M(i) is the mass of the i^thgravitating object or planet, r(i) is the x, y, and z components of the distance between the spacecraft and the i^thgravitating object or planet, g is the acceleration of gravity at the Earth's surface, and I_spis the specific impulse of the spacecraft's engine. The product gI_spis the exhaust velocity of the rocket engine.

[0133] The state equation (3) allows for the possibility that the initial (or alternately the terminal) state x(t=t₀) may be a function of the static parameter vector w. FIG. 4 provides an example when the initial state is a function of w. In this example, w[1] is defined to be the launch date from Earth. Thus, the initial spacecraft position and velocity is a function of where the Earth is in it's orbit around the Sun at the launch date w[1].

[0134] The user selected function T is not required to be continuous in time. This freedom allows the exact representation of different engine stages. An example is when an electric propulsion engine stage follows a chemical propulsion stage on a single spacecraft. In the example indicated by Equation (12), different engine stages could be represented by making I_spa discontinuous function of time. Each engine type has a characteristic value of I_sp.

[0135] In addition to providing the useful objective function 413, the state equation 414, and the constraints 415; the user of the invention may elect to include dynamic limitations 416 in the form of Equation (5). An example of a dynamic limitation is the time spacing between which navigational commands can be sent to control a spacecraft. For example, it may be possible to communicate with a spacecraft only once every ten days. Thus, the spacecraft thrust can only be changed once every ten days. For the example indicated by FIG. 4, the dynamic limitation 416 can be written, $\begin{array}{cc}v\left(t\right)=\left|\begin{array}{cc}{u}_1& \mathrm{for}t=t_0\mathrm{to}{t}_1\\ u_2& \mathrm{for}t=t_1\mathrm{to}{t}_2\\ ⋮& ⋮\\ u_N& \mathrm{for}t=t_{N-1}\mathrm{to}{t}_N,\end{array}\right|& \left(13\right)\end{array}$

[0136] where the period functions 𝑓(u_i, w, t)=u_i.

[0137] Once a useful objective function, a state equation, and any constraints and/or dynamic limitations are provided. the SDC algorithm 404 generates the optimal spacecraft trajectory, spacecraft navigation command sequence, and/or spacecraft design parameters. For the example indicated by FIG. 4, the output 417 of the SDC process is an optimal launch date w[1], an optimal flight time w[2], and an optimal thrust sequence v(t) to reach a target such that the fuel use is minimized and the launch date is after the new year of 2001.

[0138] To provide further understanding of the present invention, two numerical examples of the SDC process applied to a simple spacecraft trajectory optimization problem are given in this subsection. The numerical examples are refinements of the example outlined in FIG. 4 discussed supra. Relatively simple examples are presented for clarity. Application of the SDC process is not limited to simple engineering problems.

[0139] The goal in both numerical examples is to generate optimal spacecraft launch dates, flight times, and thrust sequences to arrive at a target near the orbit of the planet Saturn. The objective is to minimize the amount of fuel required for the trip. The optimal launch date, flight time, and thrust sequence are defined to be those which require the least fuel. In the first example, the spacecraft is assumed to have a nuclear electric propulsion (NEP) engine. NEP engines are ion propulsion engines which operate on electricity produced by a nuclear fission reactor. In the second example, the spacecraft is assumed to have a solar electric propulsion (SEP) engine. SEP engines are ion propulsion engines which operate on electricity produced by solar cell arrays. Unlike NEP engines, SEP engines can only operate near the sun where solar cell arrays are efficient. SEP engines typically only operate out to about 25% of the distance from the sun to Saturn. SEP spacecraft must coast and additionally use the gravitational fields of planets during flybys to be accelerated to higher velocities to reach outer solar system destinations like Saturn.

[0140] In both numerical examples, the state vector x(t) is defined by Equation (10), the dynamic control vector v(t) by Equation (11), and the state equation by Equation (12). The x-axis is defined to be the line between the Sun and Earth at the instant of the launch. The y-axis lies in the ecliptic plane. The z-axis points in the direction of the solar system north pole completing a right handed rectangular coordinate system. The net gravitational effect of the sun and all nine planets in the solar system is accounted for. The spacecraft is assumed to have a single ion propulsion engine with a specific impulse (measure of engine efficiency) of I_sp=3300 seconds. The static parameter vector w is defined as, $\begin{array}{cc}w=\left[\begin{array}{c}w\left[1\right]\\ w\left[2\right]\end{array}\right]=\left[\begin{array}{c}\mathrm{launch}\mathrm{date}\left(\mathrm{year}\right)\\ \mathrm{total}\mathrm{flight}\mathrm{time}\left(\mathrm{year}\right)\end{array}\right].& \left(14\right)\end{array}$

[0141] The “launch” or start of the space flight will occur from a near Earth solar orbit. The useful objective function defined by Equation (7) is used in both examples. The seventh component T[7] of the state equation (12) is the rate fuel is burned. Thus $T\left[7\right]=\frac{d{M}_s\left(t\right)}{dt},$

[0142] and the useful objective function (7) becomes $\begin{array}{cc}J={\int }_{w\left[1\right]}^{w\left[1\right]+w\left[2\right]}\frac{\Vert v\Vert }{g{I}_{\mathrm{sp}}}dt.& \left(15\right)\end{array}$

[0143] This form of the useful objective function explicitly indicates the dependence of F on the thrust vector v.

[0144] The constraints

[0145] K=x(w[1]+w[2])[1:3]−r_Saturn(w[1]+w[2])=0  (16)

[0146] L=w[1]−2001≧0  (17)

[0147] are used to require the spacecraft to reach a point 1 astronomical unit past Saturn by the end of the trajectory and to require the launch date to be after midnight Dec. 31, 2000.

[0148] In addition, the thrust magnitude is dynamically limited to be steady, allowing changes at N=100 evenly spaced intervals during the space flight. This dynamic limitation avoids calculated optimal thrust sequences v(t) which change continuously in time. This limitation reflects practical and mechanical constraints often present in spacecraft engines and control systems. Specifically $\begin{array}{cc}v\left(t\right)=\left|\begin{array}{cc}{u}_1& \mathrm{for}t=w\left[1\right]\mathrm{to}{t}_1\\ u_2& \mathrm{for}t=t_1\mathrm{to}{t}_2\\ ⋮& ⋮\\ u_N& \mathrm{for}t=t_{N-1}\mathrm{to}w\left[1\right]+w\left[2\right],\end{array}\right|& \left(18\right)\end{array}$

[0149] where the dynamic parameter vectors u_iεR³are optimized by the period formulation of the SDC process. No limitation or constraints are placed on the spacecraft trajectory that require the spacecraft to use any preprescribed flyby sequence to get to Saturn. The SDC process is allowed to select any flyby sequence.

[0150] Equations (15), (12), (16)-(17), and (18) provide the useful objective function 413, the state equation 414, the constraints 415, and the dynamic limitations 416 respectively for the NEP example. The period formulation of the Static/dynamic algorithm 404 is used to generate the optimal NEP launch date, flight time, and spacecraft trajectory.

[0151] The first step 406 in FIG. 2 of the Static/dynamic algorithm 404 incorporates the constraint equations (16)-(17) into the useful objective function J and incorporates the dynamic limitation equation (18) into the period functions 𝑓. The dynamic limitation equation (18) can be incorporated in period functions 𝑓 by defining,

[0152] 𝑓(u_i, w, t, i)=u_ifor i=1, 2, . . . , N.  (19)

[0153] The constraint equations (16) and (17) can be incorporated into the objective J using twice differentiable penalty functions. For example, the function G in the general objective equation (2) can be defined, $\begin{array}{cc}G\left(x\left(t_N\right),u_N,w,t_N,N\right)=\left\{\begin{array}{cc}{c}_l{L}^t·L& \mathrm{when}L<0\\ 0& \mathrm{when}L\ge 0\end{array}\right\}+c_d{K}^t·K& \left(20\right)\end{array}$  G(x(t_N), u_N, w, t_N, i)=0 i≠N  (21)

[0154] where c_land c_dare scaler, positive penalty function weights. The penalty weights c_land c_dare increased and the problem is resolved until no significant change in the optimal trajectory occurs with increasing values for c_land c_d. Note that any standard method for handling constraints can be used. The penalty function method is probably the simplest, and is often adequate. The useful objective with the constraints incorporated is $\begin{array}{cc}J={\int }_{w\left[1\right]}^{w\left[1\right]+w\left[2\right]}\frac{\Vert v\Vert }{g{I}_{\mathrm{sp}}}dt.+\left\{\begin{array}{cc}{c}_l{L}^t·L& \mathrm{when}L<0\\ 0& \mathrm{when}L\ge 0\end{array}\right\}+c_d{K}^t·K& \left(22\right)\end{array}$

[0155] The second step 407 in the SDC algorithm 404 in FIG. 2 is to provide arbitrary initial approximations for v(t) and w. The initial approximations of v(t) and w are used with the state equation (12) to compute x(t). Once x(t) is obtained, the useful objective value of J is calculated using the objective equation (22). To demonstrate the robustness of the SDC algorithm, the initial approximation of v(t) was set to zero for all times, $\begin{array}{cc}v\left(t\right)=\left\{\begin{array}{c}0\\ 0\\ 0\end{array}\right\}\mathrm{for}t=w\left[1\right]\mathrm{to}w\left[1\right]+w\left[2\right].& \left(23\right)\end{array}$

[0156] The initial iteration values for the launch date w[1] and flight time w[2] were arbitrarily chosen to be 12:00 am Jan. 1, 2014 and 3.25 years respectively. The initial values for v(t) and w result in the spacecraft never leaving the Earth's vicinity in its three year “flight” (the spacecraft never reaches Saturn!) The thrust sequence, Equation (23), is obviously quite poor. This example demonstrates that the SDC process is highly robust because it can be started with very poor first iteration values for v(t) and w. FIG. 5 is a plot of the first iteration spacecraft trajectory resulting from the first iteration values of v(t) and w.

[0157] The remainder of step 407 and steps 408 through 412 of the SDC algorithm 404 in FIG. 2 were implemented using SDC software on a desktop workstation. The precise algorithmic steps followed by the computer are given in subsection 5 entitled “The Steps of the SDC Algorithm For the Period Formulation.” The normalized time substitution was used which is described in detail in subsection 6 entitled “The Normalized Time Substitution.” The required number of iterations of the loop terminating in step 412 of FIG. 2 was 1589.

[0158] FIG. 6 is a plot of the optimal NEP trajectory obtained after convergence. FIG. 6 demonstrates a novel aspect of the present invention: the SDC process can identify and optimize spacecraft flybys without having the flybys prespecified. No other method for optimizing spacecraft trajectories is presently capable of doing this. As shown in FIG. 6, the SDC process identified and optimized two flybys. The first is a flyby of the planet Jupiter, which imparts kinetic energy to the spacecraft to accelerate the spacecraft toward the orbit of Saturn, thereby saving fuel. The second flyby is of the planet Saturn. The converged closest approach to Jupiter is 2.632×10⁶km. The converged closest approach to Saturn is 1.463×10⁶km. The optimal launch date (w[1]) is at 17:47:24 on Jan. 23 of the year 2012. The optimal flight time (w[2]) is 10.296 years.

[0159] FIG. 7 is a plot of the optimal thrust sequence and the velocity magnitude corresponding to the trajectory plotted in FIG. 6. The plot of thrust is given as a bar plot to emphasize that the thrust is fixed during each of the 100 periods conforming to the dynamic limitation equation (18). The thrust in the z direction does not rise substantially above zero milli-newtons on the scale of the thrust plot in FIG. 7. The magnitude of the vector sum of the x, y, and z thrusts is plotted as a staircase like line in the upper half of the thrust graph. The optimal payload fraction for this example is 64.213%. Thus 35.787% of the initial spacecraft mass must be fuel. The payload result can be scaled to any size spacecraft. The velocity magnitude plot in FIG. 7 shows the velocity gains from the Jupiter flyby around the end of year 2018 and from the Saturn flyby near year 2022. These velocity gains result in substantially less fuel being required for the trip. The calculations for Example 1, comprising 1589 iterations, were performed in 113.5 minutes on a SUN Ultra workstation.

[0160] For simplicity, Example 2 will consider the same basic problem as for the NEP case of Example 1, except that SEP engines are used instead of EP engines. Unlike NEP engines, SEP engines cannot operate beyond about 2.6 astronomical units from the sun. The target is at about 10.4 astronomical units. In order to fly SEP spacecraft beyond 2.6 astronomical units, enough kinetic energy has to be transferred to the spacecraft inside of the 2.6 astronomical unit limit. In addition, flyby assists are necessary to increase speed outside the 2.6 astronomical unit limit. This problem is significantly more difficult to solve than the NEP problem in Example 1.

[0161] The SEP engine cut-off distance introduces a new constraint on the thrust, which can be written

[0162] Trmax(r*)−∥v∥≧0,  (24)

[0163] where r* is the current distance between the spacecraft and the sun and Trmax(r*) is the maximum possible thrust which can be attained using solar energy at distance r*. A plot of the shape of Trmax(r*) used in this example is given in FIG. 8. Notice that the maximum allowed thrust Trmax(r*) drops to zero after 2.6 astronomical units.

[0164] In this example, the vector L of Equation (4) has two components $\begin{array}{cc}L=\left\{\begin{array}{c}w\left[1\right]-2001\\ \mathrm{Trmax}\left(r_{*}\right)-\Vert v\Vert \end{array}\right\}\ge \left\{\begin{array}{c}0\\ 0\end{array}\right\}.& \left(25\right)\end{array}$

[0165] The definition of L in Equation (25) will constrain the launch date w[1] and the thrust v(t) simultaneously.

[0166] Equations (15); (12); (16) and (25); and (18) provide the useful objective function 413, the state equation 414, the constraints 415, and dynamic limitations 416 respectively for the SEP example (see FIG. 4). The Static/dynamic algorithm 404 was used to generate the optimal SEP spacecraft trajectory.

[0167] The first step 406 of the Static/dynamic algorithm 404 incorporates the Equations (16) and (25) into the useful objective function J and incorporates Equation (18) into period functions 𝑓. The same approach used for the NEP case of example 1 is used here. The result is given by Equations (19), (20), and (22).

[0168] The second step 407 in the SDC algorithm 404 of FIG. 4 provides arbitrary initial approximations for v(t) and w. Then the initial approximations of v)(t) and w, and the state equation (12) are used to compute x(t). Once x(t) is obtained, the useful objective value J is calculated using the objective equation (22). To again demonstrate the robustness of the SDC algorithm, the initial value of v(t) was set to zero for all times, $\begin{array}{cc}v\left(t\right)=\left\{\begin{array}{c}0\\ 0\\ 0\end{array}\right\}\mathrm{for}t=w\left[1\right]\mathrm{to}w\left[1\right]+w\left[2\right].& \left(26\right)\end{array}$

[0169] The initial values for the launch date w[1] and flight time w[2] were arbitrarily chosen to be 12:00 am Jan. 1, 2014 and 3.25 years respectively. The initial values for v(t) and w result in the spacecraft never leaving the Earth's vicinity in its three year “flight” because the engine is never turned on. FIG. 5 is plot of the first iteration spacecraft trajectory resulting from the initial values of v(t) and w.

[0170] The remainder of step 407 and steps 408 through 412 of the SDC algorithm 404 were implemented using SDC software on a desktop workstation. The precise algorithmic steps followed by the computer are given in subsection 5 entitled “The Steps of the SDC Algorithm For the Period Formulation.” The normalized time substitution was used which is described in detail in subsection 6 entitled “The Normalized Time Substitution.” The required number of iterations of the loop terminating in step 412 of FIG. 2 was 1941.

[0171] FIG. 9 is a plot of the optimal SEP trajectory obtained after convergence. FIG. 9 demonstrates a novel aspect of the present invention: the SDC process can identify and optimize spacecraft flybys without having the flybys prespecified even when an SEP thrust constraint is in effect. No other method for optimizing spacecraft trajectories is presently capable of doing this. As shown in FIG. 9, the SDC process identified and optimized two flybys. The first is a flyby of the planet Jupiter, which imparts kinetic energy to the spacecraft to accelerate the spacecraft toward the orbit of Saturn, thereby saving fuel. The second flyby is of the planet Saturn which results in a slow down to reach the final target 1 astronomical unit outside the orbit of Saturn. The converged closest approach to Jupiter is 1.717×10⁶km. The converged closest approach to Saturn is 1.153×10⁶km. The optimal launch date (w[1]) is at 05:43:08 on Apr. 16 of the year 2012. The optimal flight time (w[2]) is 10.190 years.

[0172] FIG. 10 is a plot of the optimal thrust sequence and the velocity magnitude corresponding to the trajectory plotted in FIG. 9 for a 1040 kg initial weight spacecraft. The thrust in the z direction does not rise substantially above zero milli-newtons on the scale of the thrust plot in FIG. 10. Note that the engine is essentially shut down around year 2017 when the spacecraft exceeds the 2.6 astronomical unit distance. The optimal payload fraction for this example is 62.792%. Thus 37.208% of the initial spacecraft mass must be fuel. The velocity magnitude plot in FIG. 10 shows the velocity changes from the Jupiter flyby around the middle of year 2018 and the Saturn flyby near year 2022.

[0173] Note that despite the difficulty of sending an SEP spacecraft to a target beyond Saturn, the SDC process identified a trajectory that is similar to, and nearly as efficient, as the NEP trajectory in FIG. 6. In addition, the SDC process was started with a zero thrust initial trajectory. This demonstrates the power of the SDC process. The calculations for Example 2, comprising 1941 iterations, were performed in 141.3 minutes on a SUN Ultra workstation.

[0174] It should be understood that the present invention can be applied to virtually any spacecraft design and trajectory problem, which can differ substantially from the examples given in this subsection. While the examples given in this subsection may require certain variables, objectives and constraints, to be considered, other situations may require other variables, objectives, and constraints, in addition to, and/or in place of, those that are discussed. A detailed description of the static/dynamic control (SDC) algorithm for the period formulation is provided below. The fully continuous form of the SDC algorithm is provided in Appendix B.

[0175] The period formulation of the SDC algorithm is designed to solve problems of the following general form.

[0176] N=number of periods

[0177] n=dimension of the state vector x(t)

[0178] m=dimension of the dynamic parameter Vectors u_i(i=0, 1, . . . , N)

[0179] m′=dimension of the dynamic control vector v(t)

[0180] l=dimension of the static parameter vector w

[0181] t_i=(i=0, 1, . . . , N) times marking divisions between periods

[0182] Given a user specified useful objective function of the form, $\begin{array}{cc}J\left(u_1,u_2,\dots ,u_N,w\right)={\int }_{t_0}^{t_N}F\left(x\left(t\right),v\left(t\right),w,t\right)dt+\sum _{i=1}^{N}G\left(x\left(t_i\right),u_i,w,t_i,i\right),& \left(27\right)\end{array}$

[0183] where either of t₀<t_Nand t₀>t_Nis permitted. The period formulation of the present invention finds dynamic parameter vectors u₁, u₂, . . . , u_N, and static parameter vector w such that J is minimized (or analogously maximized),

[0184] J*=min_u,wJ(u₁, u₂, . . . , u_N, w).  (28)

[0185] The minimization (or maximization) is subject to a user specified state equation: $\begin{array}{cc}\frac{dx\left(t\right)}{dt}=T\left(x\left(t\right),v\left(t\right),w,t\right)x\left(t=t_0\right)=\Gamma \left(w\right).& \left(29\right)\end{array}$

[0186] The dynamic control vector, v(t)εR^m′, is dynamically limited to be a parametric function of time over each of N specified time intervals called “periods:” $\begin{array}{cc}v\left(t\right)=\left\{\begin{array}{cc}f\left(u_1,w,t,1\right)& \mathrm{for}t=t_0\mathrm{to}{t}_1\\ f\left(u_2,w,t,2\right)& \mathrm{for}t=t_1\mathrm{to}{t}_2\\ ⋮& ⋮\\ f\left(u_N,w,t,N\right)& \mathrm{for}t=t_{N-1}\mathrm{to}{t}_N.\end{array}\right\}& \left(30\right)\end{array}$

[0187] The dynamic control vector v(t) is a piece-wise continuous vector function of time defining the dynamic control values at each time instant. The value of v(t) may change discontinuously at period interfaces t₁, t₂, t₃, . . . , t_N−1. The functions 𝑓(u_i, w, t, i)εR^m′, i=1, 2, . . . , N are parameterized by u_iεR^mand wεR^l. The “dynamic parameter vector” for the i^thperiod is defined to be the vector of parameters u_i. The time intervals (t_i−1, t_i), i=1, 2, . . . , N and the number of periods N may be chosen arbitrarily. The elements of the static parameter vector w are fixed-in-time decision variables (i.e. non-dynamic). This formulation requires that both classes of variables, time changing (dynamic) variables and fixed-in-time (static) variables, are optimized simultaneously. The system state at time t₀, Γ(w)εRⁿis assumed to be a continuously differentiable function of the static parameter vector w.

[0188] The state vector x(t)εRⁿis a continuous function of time. It is assumed that the user selected functions 𝑓εR^m′; F and GεR; and TεRⁿare continuously differentiable with respect to state vector x and the parameter vectors u_iand w. The functions 𝑓, F, and T may be discontinuous functions of time t.

[0189] In addition, the user may define linear or nonlinear constraints of the form:

[0190] L(x(t), v(t), w, t),≧0 and/or K(x(t), v(t), w, t)=0.  (31)

[0191] Equation (31) constitute constraints relating x(t), v(t), w, and t. Any such constraint may be in the form of a constraint on a subset of x(t), v(t), w, and t. An example of such a constraint is a relationship between v(t) and t, such as v(t)≦At where A is a specified constant. Another example is a condition on v(t) alone, such as v(t)<v_maxwhere v_maxis a specified constant.

[0192] This subsection outlines the steps for implementing the SDC algorithm for the period formulation. The theoretical basis for the period formulation is presented in subsection 7.

[0193] With reference to FIG. 2, Step 406 incorporates the constraint equations (31) into the objective function (27) using a twice differentiable penalty function or other standard constraint method. Next, the dynamic limitations, if present, are incorporated into period functions 𝑓.

[0194] Step 407 provides an initial (arbitrary) value for the dynamic parameter vectors u₁, u₂, . . . , u_N, and the static parameter vector w. The initial values are denoted as {overscore (u)}₁, {overscore (u)}₂, . . . , {overscore (u)}_N, and {overscore (w)}. Next, the time evolution of the state vector x(t), subject to x(t=t₀)=Γ(w), is calculated using {overscore (u)}₁, {overscore (u)}₂, . . . , {overscore (u)}_N, and {overscore (w)} by integrating the state equation (29). The resulting state vector is denoted {overscore (x)}(t). The useful objective value J({overscore (u)}₁, {overscore (u)}₂, . . . , {overscore (u)}_N, {overscore (w)}) is calculated using Equation (27). The resulting objective value is denoted {overscore (J)}.

[0195] Step 408 initializes the program parameters P=0^n×n, Q=0^n×1, ΔR=0, S=0^l×1, Y=0^n×l, W=0^l×l, and ε=1.

[0196] Step 409 completes the following calculation for i=N, N−1, N−2, . . . , 1 comprising a first substep followed by a second substep. Since period N encompasses times between t_N−1and t_N, and since period 1 encompasses times between t₀an t₁, the calculations of step 409 proceed backward in time if t₀<t_Nand forward in time if t₀>t_N.

[0197] The first substep of step 409 is the integration of the system of coupled ordinary differential equations (32)-(40) from t=t_ito t=t_i−1using the condition equations (41)-(49). The integration results in J_x(t_i−1), J_u(t_i−1), J_w(t_i−1), J_xx(t_i−1), J_xu(t_i−1), J_xw(t_i−1), J_uu(t_i−1), J_wu(t_i−1), and J_ww(t_i−1).

[0198] {overscore ({dot over (J)})}_x=−F_x−T_x^t{overscore (J)}_x  (32)

[0199] {overscore ({dot over (J)})}_u=−F_u−T_u^t{overscore (J)}_x  (33)

[0200] {overscore ({dot over (J)})}_w=−F_w−T_w^t{overscore (J)}_x  (34)

[0201] $\begin{array}{cc}{\stackrel{\stackrel{.}{\_}}{J}}_{\mathrm{xx}}=-F_{\mathrm{xx}}-{\stackrel{\_}{J}}_{\mathrm{xx}}{T}_x-T_x^t{\stackrel{\_}{J}}_{\mathrm{xx}}-\sum _{k=1}^n{\stackrel{\_}{J}}_x\left[k\right]T_{\mathrm{xx}}\left[k,:,:\right]& \left(35\right)\\ {\stackrel{\stackrel{.}{\_}}{J}}_{\mathrm{xu}}=-F_{\mathrm{xu}}-{\stackrel{\_}{J}}_{\mathrm{xx}}{T}_u-T_x^t{\stackrel{\_}{J}}_{\mathrm{xu}}-\sum _{k=1}^n{\stackrel{\_}{J}}_x\left[k\right]T_{\mathrm{xu}}\left[k,:,:\right]& \left(36\right)\\ {\stackrel{\stackrel{.}{\_}}{J}}_{\mathrm{xw}}=-F_{\mathrm{xw}}-{\stackrel{\_}{J}}_{\mathrm{xx}}{T}_w-T_x^t{\stackrel{\_}{J}}_{\mathrm{xw}}-\sum _{k=1}^n{\stackrel{\_}{J}}_x\left[k\right]T_{\mathrm{xw}}\left[k,:,:\right]& \left(37\right)\\ {\stackrel{\stackrel{.}{\_}}{J}}_{\mathrm{uu}}=-F_{\mathrm{uu}}-{\stackrel{\_}{J}}_{\mathrm{xu}}^t{T}_u-T_u^t{\stackrel{\_}{J}}_{\mathrm{xu}}-\sum _{k=1}^n{\stackrel{\_}{J}}_x\left[k\right]T_{\mathrm{uu}}\left[k,:,:\right]& \left(38\right)\\ {\stackrel{\stackrel{.}{\_}}{J}}_{\mathrm{wu}}=-F_{\mathrm{wu}}-{\stackrel{\_}{J}}_{\mathrm{xw}}^t{T}_u-T_w^t{\stackrel{\_}{J}}_{\mathrm{xu}}-\sum _{k=1}^n{\stackrel{\_}{J}}_x\left[k\right]T_{\mathrm{wu}}\left[k,:,:\right]& \left(39\right)\\ {\stackrel{\stackrel{.}{\_}}{J}}_{\mathrm{ww}}=-F_{\mathrm{ww}}-{\stackrel{\_}{J}}_{\mathrm{xw}}^t{T}_w-T_w^t{\stackrel{\_}{J}}_{\mathrm{xw}}-\sum _{k=1}^n{\stackrel{\_}{J}}_x\left[k\right]T_{\mathrm{ww}}\left[k,:,:\right]& \left(40\right)\end{array}$

[0202] Equations (32)-(40) are respectively subject to the conditions at time t_i:

[0203] {overscore (J)}_x(t_i)=G_x({overscore (x)}(t_i), {overscore (u)}_i, {overscore (w)}, t_i)+Q  (41)

[0204] {overscore (J)}_u(t_i)=G_u({overscore (x)}(t_i), {overscore (u)}_i, {overscore (w)}, t_i)  (42)

[0205] {overscore (J)}_w(t_i)=G_w({overscore (x)}(t_i), {overscore (u)}_i, {overscore (w)}, t_i)+S  (43)

[0206] {overscore (J)}_xx(t_i)=G_xx({overscore (x)}(t_i), {overscore (u)}_i, {overscore (w)}, t_i)+P  (44)

[0207] {overscore (J)}_xu(t_i)=G_xu({overscore (x)}(t_i), {overscore (u)}_i, {overscore (w)}, t_i)  (45)

[0208] {overscore (J)}_xw(t_i)=G_xw({overscore (x)}(t_i), {overscore (u)}_i, {overscore (w)}, t_i)+Y  (46)

[0209] {overscore (J)}_uu(t_i)=G_uu({overscore (x)}(t_i), {overscore (u)}_i, {overscore (w)}, t_i)  (47)

[0210] {overscore (J)}_wu(t_i)=G_wu({overscore (x)}(t_i), {overscore (u)}_i, {overscore (w)}, t_i)  (48)

[0211] {overscore (J)}_ww(t_i)=G_ww({overscore (x)}(t_i), {overscore (u)}_i, {overscore (w)}, t_i)+W  (49)

[0212] where {overscore (J)}_x(t)εRⁿ, {overscore (J)}_u(t)εR^m, {overscore (J)}_w(t)εR^l, {overscore (J)}_xx(t)εR^n×n, {overscore (J)}_xu(t)εR^n×m, {overscore (J)}_xw(t)εR^n×l, {overscore (J)}_wu(t)εR^l×m, {overscore (J)}_ww(t)εR^l×l, and {overscore (J)}_uu(t)εR^m×m. The superscript t indicates the transpose operation. All derivatives of F, G, and T in Equations (32)-(49) are evaluated using {overscore (x)}(t), {overscore (u)}_i, and {overscore (w)} and can be evaluated either analytically or numerically. The condition equations (41), (42), . . . , and (49) apply to Equations (32), (33), . . . , and (40) respectively.

[0213] The second substep of step 409 is implemented as follows. If i≧2, then the results of the integration of Equations (32)-(40) are used to evaluate and store:

[0214] α_i=−{overscore (J)}_uu⁻¹(t_i−1){overscore (J)}_u(t_i−1)

[0215] β_i=−{overscore (J)}_uu⁻¹(t_i−1){overscore (J)}_xu^t(t_i−1)

[0216] γ_i=−{overscore (J)}_uu⁻¹(t_i−1){overscore (J)}_wu^t(t_i−1).

[0217] Additionally, the values of ΔR, P, Q, S, W, and Y are replaced according to:

[0218] ΔR=ΔR+{overscore (J)}_u^t(t_i−1){overscore (J)}_uu⁻¹(t_i−1){overscore (J)}_u(t_i−1)=ΔR−{overscore (J)}_u^t(t_i−1)α_i

[0219] P={overscore (J)}_xx(t_i−1)−{overscore (J)}_xu(t_i−1){overscore (J)}_uu⁻¹(t_i−1){overscore (J)}_xu^t(t_i−1)={overscore (J)}_xx(t_i−1)+{overscore (J)}_xu(t_i−1)β_i

[0220] Q={overscore (J)}_x(t_i−1)−{overscore (J)}_xu(t_i−1){overscore (J)}_uu⁻¹(t_i−1){overscore (J)}_u(t_i−1)={overscore (J)}_x(t_i−1)+β_i^t{overscore (J)}_u(t_i−1)

[0221] S={overscore (J)}_w(t_i−1)−{overscore (J)}_wu(t_i−1){overscore (J)}_uu⁻¹(t_i−1){overscore (J)}_u(t_i−1)={overscore (J)}_w(t_i−1)+γ_i^t{overscore (J)}_u(t_i−1)

[0222] W={overscore (J)}_ww(t_i−1)−{overscore (J)}_wu(t_i−1){overscore (J)}_uu⁻¹(t_i−1){overscore (J)}_wu^t(t_i−1)={overscore (J)}_ww(t_i−1)+{overscore (J)}_wu(t_i−1)γ_i

[0223] Y={overscore (J)}_xw(t_i−1)−{overscore (J)}_xu⁻¹(t_i−1){overscore (J)}_uu⁻¹(t_i−1){overscore (J)}_wu^t(t_i−1)={overscore (J)}_xw(t_i−1)+β_i^t{overscore (J)}_wu^t(t_i−1)

[0224] If i=1 then the results of the integration of Equations (32)-(40) are used to evaluate and store: ${\alpha }_1={\left({\stackrel{\_}{J}}_{\mathrm{uu}}\left(t_0\right)-{\stackrel{\widetilde{\_}}{J}}_{\mathrm{wu}}^t\left(t_0\right){\stackrel{\widetilde{\_}}{J}}_{\mathrm{ww}}^{-1}\left(t_0\right){\stackrel{\widetilde{\_}}{J}}_{\mathrm{wu}}\left(t_0\right)\right)}^{-1}\left({\stackrel{\widetilde{\_}}{J}}_{\mathrm{wu}}^t\left(t_0\right){\stackrel{\widetilde{\_}}{J}}_{\mathrm{ww}}^{-1}\left(t_0\right){\stackrel{\widetilde{\_}}{J}}_w\left(t_0\right)-{\stackrel{\widetilde{\_}}{J}}_u\left(t_0\right)\right),\xi =-{\stackrel{\widetilde{\_}}{J}}_{\mathrm{ww}}^{-1}\left(t_0\right)\left({\stackrel{\widetilde{\_}}{J}}_w\left(t_0\right)+{\stackrel{\_}{J}}_{\mathrm{wu}}\left(t_0\right){\alpha }_1\right).$

[0225] and the value of ΔR is replaced according to: $\Delta R=\Delta R+{\stackrel{\widetilde{\_}}{J}}_w^t\left(t_0\right){\stackrel{\widetilde{\_}}{J}}_{\mathrm{ww}}^{-1}\left(t_0\right){\stackrel{\widetilde{\_}}{J}}_w\left(t_0\right)+\left({\stackrel{\widetilde{\_}}{J}}_w^t\left(t_0\right){\stackrel{\widetilde{\_}}{J}}_{\mathrm{ww}}^{-1}\left(t_0\right){\stackrel{\widetilde{\_}}{J}}_{\mathrm{wu}}\left(t_0\right)-{\stackrel{\_}{J}}_u^t\left(t_0\right)\right){\alpha }_1$

[0226] where ${\stackrel{\widetilde{\_}}{J}}_w^t\stackrel{.}{=}{\stackrel{\_}{J}}_x^t\left(t_0\right){\stackrel{\_}{\Gamma }}_w+{\stackrel{\_}{J}}_w^t\left(t_0\right)$ ${\stackrel{\widetilde{\_}}{J}}_{\mathrm{ww}}\stackrel{.}{=}2{\stackrel{\_}{J}}_x^t\left(t_0\right){\stackrel{\_}{\Gamma }}_{\mathrm{ww}}+{\stackrel{\_}{\Gamma }}_w^t{\stackrel{\_}{J}}_{\mathrm{xx}}\left(t_0\right){\stackrel{\_}{\Gamma }}_w+{\stackrel{\_}{J}}_{\mathrm{ww}}\left(t_0\right)+{\stackrel{\_}{\Gamma }}_w^t{\stackrel{\_}{J}}_{\mathrm{xw}}\left(t_0\right)+{\stackrel{\_}{J}}_{\mathrm{xw}}^t\left(t_0\right){\stackrel{\_}{\Gamma }}_w$ ${\stackrel{\widetilde{\_}}{J}}_{\mathrm{wu}}\stackrel{.}{=}{\stackrel{\_}{\Gamma }}_w^t{\stackrel{\_}{J}}_{\mathrm{xu}}\left(t_0\right)+{\stackrel{\_}{J}}_{\mathrm{wu}}\left(t_0\right).$

[0227] Step 410 is the calculation of the updated static and dynamic parameter vectors, w and v(t), and the corresponding state vector x(t) for each time period in the time order of time periods 1, 2, . . . , j, . . . , N as follows: For j=1, the updated dynamic and static parameter vectors, u₁and w, are respectively calculated via:

[0228] u₁={overscore (u)}₁+εα₁  (50)

[0229] w={overscore (w)}+εξ.  (51)

[0230] The resulting values u₁, w and the updated state x(t₀)=Γ(w) are used to integrate the state equation (29) from t₀to t₁. For j=2, 3, . . . , N the updated dynamic parameter vector u_jis calculated using:

[0231] u_j={overscore (u)}_j+εα_j+β_j(x(t_j−1)−{overscore (x)}(t_j−1))+εγ_jξ  (52)

[0232] The resulting value of u_jis used to integrate the state equation (29) from t_j−1to t_j.

[0233] Concurrently with Step 410, the useful objective function (27) is evaluated using the updated control. The resulting objective value is denoted J(ε).

[0234] Step 411 is a comparison of the useful objective value J(ε) with the useful objective value {overscore (J)}. If there is satisfactory change, the values of the variables {overscore (J)}, {overscore (u)}, {overscore (w)}, and {overscore (x)}(t) are respectively replaced by the values of J(ε), u, w, and x(t). A satisfactory change is a satisfactory reduction to J(ε) if J is to be minimized. A satisfactory change is a satisfactory increase to J(ε) if J is to be maximized. If there is not a satisfactory change, then ε is reduced and Step 410 is repeated. For example, ε can be reduced by one half. A satisfactory reduction criterion for minimization that can be used is: $\begin{array}{cc}\stackrel{\_}{J}-J\left(\epsilon \right)\ge \frac{1}{2}\left(\frac{{\epsilon }^2}{2}-\epsilon \right)\Delta R& \left(53\right)\end{array}$

[0235] A satisfactory increase criterion for maximization that can be used is: $\begin{array}{cc}J\left(\epsilon \right)-\stackrel{\_}{J}\ge -\frac{1}{2}\left(\frac{{\epsilon }^2}{2}-\epsilon \right)\Delta R& \left(54\right)\end{array}$

[0236] The criteria (53) and (54) are simply requirements that the observed improvement in the useful objective be at least half the improvement that is expected if the problem were exactly quadratic (see subsection 8 titled “The Expected Quadratic Change for the Period Formulation”).

[0237] Step 412 is a test for convergence. If convergence is achieved, then the calculation ends. If convergence is not yet achieved, the Steps 408, 409, 410, and 411 are repeated. A convergence test that can be used is: convergence is achieved if the update law equations (50)-(52) with ε=1 result in 1) a change in the objective that is very small (i.e., within a specified tolerance) and is nearly equal to −½ΔR, and 2) {tilde over (J)}_ww(t₀); J_uu(t_i) for i=0, 2, 3, . . . , N−1; and ${\stackrel{\_}{J}}_{\mathrm{uu}}\left(t_0\right)-{\widetilde{{\stackrel{\_}{J}}^t}}_{\mathrm{wu}}\left(t_0\right){\stackrel{\_}{J}}_{\mathrm{ww}}^{-1}\left(t_0\right){\widetilde{\stackrel{\_}{J}}}_{\mathrm{wu}}\left(t_0\right)$

[0238] are all positive definite for minimization or all negative definite for maximization. This convergence test satisfies the necessary and sufficient conditions of optimality. Note, also, that when any of the above arc not positive definite (or negative definite for maximization), a variation of the well known Trust Region method or shifting method can be used.

[0239] This subsection presents a method to solve dynamic optimization problems with variable beginning and ending times and variable time horizons. The goal is to optimize the trajectory simultaneously with optimizing the time t₀, and length of time |t_N−t₀| used for the trajectory. The time t₀is the initial time if t₀<t_Nor it is the final time if t₀>t_N. The general method presented here was originally developed for the specific problem of space flight trajectory optimization but is applicable to general physical problems. The method developed here uses the strengths of the SDC algorithm to simultaneously optimize static and dynamic variables.

[0240] It is not immediately possible for SDC to optimize the start-time and the time-horizon. The SDC algorithm requires a fixed start and end to the integration. A further complication is that if the time-horizon t_Nt₀to changes, then individual periods (t_i−1, t_i) must also change in length. In particular, the relationships t_i−1<t_iis must be maintained if t_N>t₀or the relationships t_i−1>t_imust be maintained if t_N<t₀. The method in this subsection assumes that as the time-horizon t_Nto t₀changes in size, the sizes of the individual periods t_ito t_i−1change such that t_i−t_i−1/t_N−t₀remains constant. For example, if |t_N−t₀| is halved, then each period is halved. This method guarantees that the relationships t_i−1<t_ior, alternately, t_i−1>t_iare maintained. It is easy to generalize the method given in this subsection to optimize the time-duration of all periods (t_i−1, t_i) simultaneously by assigning a component of w or u_ito each period.

[0241] The normalized time substitution addresses the following general problem: $\begin{array}{cc}{J}^{*}=\underset{u_i,w}{\min }J=\underset{u_i,w}{\min }{\int }_{w\left[1\right]}^{w\left[1\right]+w\left[2\right]}F\left(x\left(t\right),v\left(t\right),w,t\right)dt+\sum _{N=i}^{1}G\left(x\left(t_i\right),u_i,w,t_i,i\right),& \left(55\right)\end{array}$

[0242] subject to the equation of state: $\begin{array}{cc}\frac{dx\left(t\right)}{dt}=T\left(x\left(t\right),v\left(t\right),w,t\right)x\left(t=w\left[1\right]\right)=\Gamma \left(w\right),& \left(56\right)\end{array}$

[0243] where v(t) is subject to the dynamic limitation: $\begin{array}{cc}v\left(t\right)=\left|\begin{array}{cc}f\left(u_1,w,t,1\right)& \mathrm{for}t=t_0\mathrm{to}{t}_1\\ f\left(u_2,w,t,2\right)& \mathrm{for}t=t_1\mathrm{to}{t}_2\\ ⋮& ⋮\\ f\left(u_N,w,t,N\right)& \mathrm{for}t=t_{N-1}\mathrm{to}{t}_N,\end{array}\right|& \left(57\right)\end{array}$

[0244] and where x, v, and w are subject to general constraints of the form:

[0245] L(x, v, w, t)≧0 and/or K(x, v, w, t)=0.  (58)

[0246] The static parameter vector w is assumed to have at least two components. The first two components are w[1] the start time of the problem and w[2] the duration of the problem. Note that Equations (55)-(57) are not in a form which SDC can solve directly. Equations (55)-(57) can be transformed into a mathematically equivalent problem which SDC can solve directly by making the substitution:

[0247] t=w[1]+w[2]τ 0≦τ≦1  (59)

[0248] where t is the “real time” or clock time and τ is the “normalized time.”

[0249] If Equation (59) is substituted into Equation (55) then: $\begin{array}{cc}\underset{u_i,w}{\min }{\int }_0^{1}F\left(x\left(w\left[1\right]+w\left[2\right]\tau \right),v\left(w\left[1\right]+w\left[2\right]\tau \right),w,w\left[1\right]+w\left[2\right]\tau \right)w\left[2\right]d\tau +\sum _{i=1}^N{G}_i\left(x\left(w\left[1\right]+w\left[2\right]{\tau }_i\right),u_i,w,w\left[1\right]+w\left[2\right]{\tau }_i,i\right),& \left(60\right)\end{array}$

[0250] where ${\tau }_i=\frac{t_i-w\left[1\right]}{w\left[2\right]}$

[0251] and dt=w[2]dτ. Depending on the physical interpretations of x, v, and w, the values x, v, and w may depend on w[2]. This may be understood by considering the example when x is equal to velocity. If the value of x is 5 in real time units or t-time, then, if w[2]=2, x is w[2]×5=10 in the normalized time units or τ-time. The variables {tilde over (x)}, {tilde over (v)}, and {tilde over (w)} are defined to be the state vector, dynamic control vector, and static parameter vector, respectively, in τ-time units. In general,

[0252] x(t=w[1]+w[2]τ)=(w[2])^nx{tilde over (x)}(τ)  (61)

[0253] v(t=w[1]+w[2]τ)=(w[2])^nv{tilde over (v)}(τ)  (62)

[0254] w=(w[2])^nw{tilde over (w)}.  (63)

[0255] For example, if x is a velocity, then nx=−1. In general, if the units of x include time^nxthen:

[0256] x(t=w[1]+w[2)τ=(w2])^nx{tilde over (x)}(τ)

[0257] because $w\left[2\right]=\frac{t\mathrm{time}\mathrm{units}}{\tau \mathrm{time}\mathrm{units}}$

[0258] is the conversion factor. If the units of x do not include time then nx=0. The same arguments hold for v and nv. For example, if v is interpreted as thrust, then the units of v are [mass][length][time]⁻². It follows that the conversion is v(t=w[1]+w[2]τ)=(w[2])⁻²{tilde over (v)}(τ). For w and nw, note that the first two components of w define the substitution and hence do not require conversion (i.e. nw=0 for w[1] and w[2].) If w has more than two components, then nw will be non-zero for components of w with units of time. Note: in general, different components of state and control vectors may have different units and hence require different factors of w[2] for the conversion. An example is when the first 3 components of x are position coordinates (units of length no time units) and the second set of three components are velocity (units of length and time⁻¹.) For this example, the conversion for the first three components is (w[2])⁰=1 and the conversion for the second three components is ${\left(w\left[2\right]\right)}^{-1}=\frac{1}{w\left[2\right]}.$

[0259] Thus (w[2])^nx, (w[2])^nv, and (w[2])^nware interpreted as vectors εRⁿ, R^m, and R^mrespectively. For example, ${\left(w\left[2\right]\right)}^{\mathrm{nx}}=\left[\begin{array}{c}{\left(w\left[2\right]\right)}^{\mathrm{nx}\left[1\right]}\\ {\left(w\left[2\right]\right)}^{\mathrm{nx}\left[2\right]}\\ ⋮\\ {\left(w\left[2\right]\right)}^{\mathrm{nx}\left[n\right]}\end{array}\right]$

[0260] where nx[i] is the appropriate factor of w[2] for the i^thcomponent of the state x. Further, products like (w[2])^nx{tilde over (x)}(τ) are interpreted as an element by element product εRⁿ: ${\left(w\left[2\right]\right)}^{\mathrm{nx}}\tilde{x}\left(\tau \right)=\left[\begin{array}{c}{\left(w\left[2\right]\right)}^{\mathrm{nx}\left[1\right]}\widetilde{x_1}\left(\tau \right)\\ {\left(w\left[2\right]\right)}^{\mathrm{nx}\left[2\right]}\widetilde{x_2}\left(\tau \right)\\ ⋮\\ {\left(w\left[2\right]\right)}^{\mathrm{nx}\left[n\right]}\widetilde{x_n}\left(\tau \right)\end{array}\right]$

[0261] Using the definitions for {tilde over (x)}, {tilde over (v)}, and {tilde over (w)}, Equation (60) can be rewritten in terms of the state and control in τ-time units: $\begin{array}{cc}\underset{u_i,w}{\min }{\int }_0^{1}F\left({\left(w\left[2\right]\right)}^{\mathrm{nx}}\tilde{x}\left(\tau \right),{\left(w\left[2\right]\right)}^{\mathrm{nv}}\tilde{v}\left(\tau \right),{\left(w\left[2\right]\right)}^{\mathrm{nw}}w,w\left[1\right]+w\left[2\right]\tau \right)w\left[2\right]d\tau +\sum _{i=1}^N{G}_i\left({\left(w\left[2\right]\right)}^{\mathrm{nx}}\tilde{x}\left({\tau }_i\right),u_i,{\left(w\left[2\right]\right)}^{\mathrm{nw}}\tilde{w},w\left[1\right]+w\left[2\right]{\tau }_i,i\right).& \left(64\right)\end{array}$

[0262] The time transformation of Equation (59) will be used to transform the state equation, the dynamic limitation, and the constraints into forms that explicitly relate to the normalized time variable τ.

[0263] Substituting Equation (59) into the continuous-time state equation (56) yields: $\begin{array}{cc}\frac{dx\left(t=w\left[1\right]+w\left[2\right]\tau \right)}{w\left[2\right]dt}=T\left(x\left(w\left[1\right]+w\left[2\right]\tau \right),v\left(w\left[1\right]+w\left[2\right]\tau \right),w,t=w\left[1\right]+w\left[2\right]\tau \right)x\left(t=w\left[1\right]\right)=\Gamma \left(w\right).& \left(65\right)\end{array}$

[0264] Using the definition equations (61)-(63) Equation (65) can be rewritten as: $\begin{array}{cc}\frac{d\left({\left(w\left[2\right]\right)}^{\mathrm{nx}}\tilde{x}\left(\tau \right)\right)}{w\left[2\right]d\tau }=T\left({\left(w\left[2\right]\right)}^{\mathrm{nx}}\tilde{x}\left(\tau \right),{\left(w\left[2\right]\right)}^{\mathrm{nv}}\tilde{v}\left(\tau \right),{\left(w\left[2\right]\right)}^{\mathrm{nw}}\tilde{w},w\left[1\right]+w\left[2\right]\tau \right),& \left(66\right)\end{array}$

[0265] subject to the condition:

[0266] (w[2])^nx{tilde over (x)}(τ=0)=Γ((w[2])^nw{tilde over (w)}).  (67)

[0267] If the factors of w[2] are shifted to the right-hand-sides of Equations (66) and (67) then: $\begin{array}{cc}\frac{d\tilde{x}\left(\tau \right)}{d\tau }={\left(w\left[2\right]\right)}^{1-\mathrm{nx}}T\left({\left(w\left[2\right]\right)}^{\mathrm{nx}}\tilde{x}\left(\tau \right),{\left(w\left[2\right]\right)}^{\mathrm{nv}}\tilde{v}\left(\tau \right),{\left(w\left[2\right]\right)}^{\mathrm{nw}}\tilde{w},w\left[1\right]+w\left[2\right]\tau \right),& \left(68\right)\end{array}$

[0268] subject to the condition:

[0269] {tilde over (x)}(τ=0)=(w[2])^−nxΓ((w[2])^nw{tilde over (w)}),

[0270] which is the state equation for {tilde over (x)}.

[0271] Substitution of Equation (59) into the dynamic limitation equation (57) yields: $\begin{array}{cc}v\left(w\left[1\right]+w\left[2\right]\tau \right)=\left|\begin{array}{cc}f\left(u_1,w,w\left[1\right]+w\left[2\right]\tau ,1\right)& \mathrm{for}\tau =0\mathrm{to}\frac{t_1-w\left[1\right]}{w\left[2\right]}\\ f\left(u_1,w,w\left[1\right]+w\left[2\right]\tau ,2\right)& \mathrm{for}\tau =\frac{t_1-w\left[1\right]}{w\left[2\right]}\mathrm{to}\frac{t_2-w\left[1\right]}{w\left[2\right]}\\ ⋮& ⋮\\ f\left(u_N,w,w\left[1\right]+w\left[2\right]\tau ,N\right)& \mathrm{for}\tau ={\mathrm{fract}}_{N-1}-w\left[1\right]w\left[2\right]\mathrm{to}1.\end{array}\right|& \left(69\right)\end{array}$

[0272] If τ_iis defined to be ${\tau }_i=\frac{t_i-w\left[1\right]}{w\left[2\right]}$