Patent Publication Number: US-7216004-B2

Title: Method and system for optimization of general problems

Description:
CROSS-REFERENCES TO RELATED APPLICATIONS 
   This application is a continuation of U.S. application Ser. No. 10/745,394, filed Dec. 23, 2003, now U.S. Pat. No. 7,072,723, which is a continuation-in-part application of U.S. application Ser. No. 10/693,729, filed Oct. 23, 2003, which claims the benefit of provisional patent application Nos. 60/420,921 and 60/420,922, both filed Oct. 23, 2002. 

   TECHNICAL FIELD 
   The present invention is related to general optimization methods employed to determine relatively optimal solutions for various types of problems and to control of computational processes and construction of computational approaches to solving complex problems and, in particular, to a general computational method and system for minimizing a general class of objective functions, for adjusting solutions to optimization over time, and for constructing and achieving relatively optimal computational control over a wide range of computational processes. 
   BACKGROUND OF THE INVENTION 
   Significant research-and-development effort is currently being applied to the mathematical and computational field of optimization. Optimization techniques are applied to a wide variety of everyday and complex theoretical problems in order to find solutions to those problems, including cost-effective and time-efficient control strategies for complex systems and organizations. For example, optimization techniques may be applied to schedule labor and equipment resources within a large manufacturing organization in order to most cost-effectively produce a given amount of manufactured product to meet a projected demand for the product at various points in time. As another example, optimization techniques may be applied to a large control problem, such as the control of traffic lights within a city, in order to allow for as large a volume of traffic to move through the city as possible without generating unreasonable delays at various intersections and for those traveling on various secondary and tertiary routes within the city. Many natural systems can be viewed as seeking relatively optimal solutions to complexly defined problems governed by physical characteristics and principles and constrained by a large number of physical constraints. For example, complex chemical reactions can be considered to be governed by rather simple, thermodynamic principles by which individual components, such as molecules, seek minimum energy states and maximum entropy, constrained by various ways in which molecules can reconfigure themselves and exchange energy amongst themselves according to quantum mechanics. Even more complex biological systems can also be considered to be governed by chemical and thermodynamic principles as well as by analogous principles involving information content and exchange, compartmentalization and modular design, and other considerations. Thus, optimization problems may encompass an extremely wide variety of mathematically and computationally expressed models of natural phenomena, design, organization, and control of complex systems, and organization, transmission, and processing of information. 
   Current approaches for finding near optimal and optimal solutions for mathematically modeled problems are limited. When the number of decision variables and constraints in such problems increases from the small number of variables and constraints normally employed in simple, textbook problems to the large number of variables and constraints normally encountered in real-world systems, the computational resources required for seeking near optimal and optimal solutions increase exponentially in most cases. Current techniques cannot be satisfactorily applied to any but the simplest types of problems. Many currently available techniques involve applying oppressive constraints to optimization models, such as requiring variables to be continuous and requiring the hyper-dimensional volume representing the set of possible solutions of optimizations problems, or problem domain, to be convex. 
   Researchers, developers, system modelers, and investigators of many different types of complex system behavior have therefore recognized the need for more generally applicable and efficient methods and systems for optimization. 
   SUMMARY OF THE INVENTION 
   Various embodiments of the present invention include methods and systems for generic optimization of problems by an approach to minimizing functions over high-dimensional domains that mathematically model problems for which near optimal or optimal solutions are sought. These embodiments receive a mathematical description of a system, in symbolic form, that includes decision variables of various types, including real-number-valued, integer-valued, and Boolean-valued decision variables. The objective function with high dimensional domain that is minimized in order to find near optimal and optimal solutions may be accompanied by a variety of constraints on the values of the decision variables, including inequality and equality constraints. Various embodiments of the present invention incorporate the objective function and constraints into a global objective function. The local and global minima of this function occur at critical points. One or more of these critical points is a near optimal or optimal solution for the modeled problem. Various embodiments of the present invention transform the global objective function and a procedure for finding critical points into a system of differential equations in terms of continuous variables and parameters, so that powerful, polynomial-time methods for solving differential equations can be applied for identifying critical points of the function. These embodiments undertake an interior-point-based method and employ a global gradient-descent field formulation in order to efficiently traverse the hyperdimensional domain of the global objective function to find local minima. Near a local critical point, where the global descent field is ineffective in providing better approximants of the variables, embodiments of the present invention provides a steering method that generates a local descent-field across the local critical point towards the solution. Once the vicinity of the local critical point is left behind, the global gradient descent method is resumed, and various embodiments of the present invention employ a method to establish a reasonable direction for resuming a descent path through the hyper-dimensional domain toward local minima. One of these local minima may represent a near optimal solution of the problem. Embodiments of the present invention also provide for distribution and decomposition of the global gradient descent-field and local gradient descent-field optimization methods using multiple threads and agents, respectively, in order to allow for parallel computation and increased time efficiency. Various embodiments of the present invention further include approaches for adjusting solutions to optimization problems relatively continuously in time, without needing to recompute the optimization solution de novo. While many embodiments of the present invention are specifically directed to various classes of optimization problems, other embodiments of the present invention provide a more general approach for constructing complex hierarchical computational processes and for optimally or near optimally controlling general computational processes. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
       FIGS. 1A–B  illustrate a very simple, two-dimensional optimization problem. 
       FIG. 2  is a control-flow diagram for a highest level view of an optimization method and system that represents one embodiment of the present invention. 
       FIG. 3  illustrates a hypothetical problem used to illustrate modeling and preprocessing. 
       FIGS. 4–7  illustrate further aspects and characteristics of the hypothetical problem using the illustration conventions of  FIG. 3 . 
       FIG. 5  illustrates costs associated with shipping goods from a manufacturing plant to warehouses and from warehouses to customers. 
       FIG. 6  illustrates additional characteristics of the hypothetical problem. 
       FIG. 7  shows a simple two-manufacturing-plant, three-potential-warehouse, five-customer hypothetical problem associated with numerical costs, capacities, and demands. 
       FIG. 8  is a more detailed, control-flow diagram for the preprocess step of the high-level control-flow diagram shown in  FIG. 2 . 
       FIG. 9  shows an initial symbolic, mathematically expressed model for the specific optimization problem illustrated in  FIG. 7 . 
       FIG. 10  shows transformation of the hypothetical mathematical model into a model employing only floating point variables. 
       FIG. 11  shows the mathematical model for the hypothetical problem transformed to standard form. 
       FIGS. 12A–D  illustrate the effect of binding variables on transformed equality constraints. 
       FIG. 13  shows the symbolic, mathematically expressed model for the hypothetical problem following transformation of equality constraints to inequality constraints. 
       FIG. 14  graphically illustrates the effect of the binding variables. 
       FIG. 15  shows the mathematical model for the hypothetical problem following introduction of the envelope variable z. 
       FIGS. 16A–D  illustrate the concept of barrier functions. 
       FIG. 17  shows the symbolic, mathematical model for the hypothetical problem following addition of barrier constraints. 
       FIG. 18  illustrates the use of a gradient field Within a problem domain in order to provide a direction for an optimization trajectory leading from an initial point to a local minimum. 
       FIG. 19  illustrates an optimization trajectory. 
       FIG. 20  graphically illustrates an iterative optimization method based on a near-optimal state-vector trajectory y* and near-optimal control policies u* and r*. 
       FIG. 21  provides additional, graphical illustration of the overall iterative optimization method. 
       FIG. 22  is a control-flow diagram that describes the iterative optimization step ( 203  in  FIG. 2 ). 
       FIG. 23  shows a surface boundary for a bowl-shaped problem domain with a number of prominent central protrusions. 
       FIG. 24  shows the same problem domain as shown in  FIG. 23  with critical points explicitly indicated. 
       FIG. 25  shows a magnified view of a critical point. 
       FIG. 26  illustrates a first-search technique. 
       FIG. 27  illustrates a second strategy for searching a problem domain. 
       FIG. 28  illustrates the modules and modular organization of an optimization thread. 
       FIG. 29  illustrates, in part, the contents of the database maintained by the super module. 
       FIGS. 30–34  illustrate, using the illustration conventions of  FIG. 28 , the flow of information among the hierarchically organized modules of an optimization thread during one outer-loop iteration of the optimization method. 
       FIG. 32  illustrates data flow during a next step of the inner loop ( 2205  in  FIG. 22 ). 
       FIG. 33  indicates that the date-exchange steps illustrated in  FIGS. 31 and 32  are successively repeated in a number of inner-loop iterations. 
       FIG. 34  illustrates the data exchange involved in costate vector calculations undertaken by the p module  2810  in step  2209  of  FIG. 22 . 
       FIG. 35  graphically illustrates computation of the w vector. 
       FIG. 36  illustrates multiplication of the Hessian matrix by the w vector to produce the G vector in a distributed fashion. 
       FIG. 37  shows a decomposed multiplication of the Hessian matrix by the w vector. 
       FIG. 38  shows a mathematical manipulation of the computation of a G-vector partition that facilitates decomposition of the Hessian-matrix inversion by y-module agents. 
       FIG. 39  illustrates computation of a w-vector partition w 1  by the second y-module agent. 
       FIG. 40  is a control-flow diagram illustrating a highest conceptual level of the super module. 
       FIG. 41  is a control-flow diagram for the routine “initial point” called in step  4004  of  FIG. 40 . 
       FIG. 42  is a control-flow diagram for the routine “choose_next_T.” 
       FIG. 43  is a control-flow diagram of the routine “evaluate,” called in step  4012  of  FIG. 40 . 
       FIG. 44  is a control-flow diagram for the routine “report.” 
       FIG. 45  is a control-flow diagram for the routine “integrate_u_and_y_to_T,” called in step  4010  of  FIG. 40 . 
       FIG. 46  is a control-flow diagram for the routine “integrate y.” 
       FIG. 47  is a control-flow diagram for the routine “compute Δτ,” called in step  4625  of  FIG. 46 . 
       FIG. 48  shows that the Hessian matrix can be decomposed by spectral decomposition. 
       FIG. 49  shows successive state-function calculations by an optimization method for a control problem. 
       FIG. 50  illustrates a shortened interval between calculations of the state functions illustrated in  FIG. 49 . 
       FIG. 51  illustrates an analogy for consequence arising from repeated, current-state-based optimizations. 
       FIG. 52  illustrates the concept of a sliding window. 
       FIG. 53  illustrates a state function calculated in a prior window extends from time T 1 −ΔT to time T T −ΔT. 
       FIG. 54  illustrates extension of the previously calculated state function to extend over the time interval of the current window. 
       FIG. 55  illustrates consideration, at a current point in time, of past events and environmental changes that may result in a previously calculated state function being out-of-date, and not optimal, at the current point in time, as well as an approach to incorporating information from the previous state function in finding a new, optimal state function for the current frame of reference represented by the current window. 
       FIG. 56  shows repair of a previously calculated state function in order to provide a near optimal state function for a current time window. 
       FIG. 57  illustrates a single level of meta-control of a computational process. 
       FIG. 58  illustrates multiple levels of meta-control and/or hierarchical construction of a complex computational process. 
   

   DETAILED DESCRIPTION OF THE INVENTION 
   The present invention relates to the generic optimization of mathematically-formulated problems, relatively continuous adjustment of solutions to optimization control problems, and to a more abstract, higher-level hierarchical construction of, an optimal or near-optimal control of computational processes. Various embodiments of the present invention include complex software programs running on a single-processor computer system, or running, in parallel, on multi-processor computer systems, or on a large number of distributed, interconnected single-and/or-multiple processor computer systems. These programs are quite mathematical in nature, and involve elaborate and difficult-to-describe mathematical techniques. For this reason, the present invention is described in part, below, with reference to a concrete problem, with reference to many hundreds of equations, and with reference to many graphical illustrations and control-flow diagrams. Although mathematical expressions, alone, may be sufficient to fully describe and characterize the programs that implement processes that represent various embodiments of the present invention to those skilled in the art of software development, the more graphical, concrete-problem-oriented, and control-flow-diagram approaches included in the following discussion are intended to illustrate various embodiments of the present invention in a variety of different ways that may be particularly accessible to readers with particular backgrounds. The present invention is described below in the following subsections: (1) Optimization; (2) Repair; and (3) Meta-Control and Hierarchical Process Construction. 
   Optimization 
     FIGS. 1A–B  illustrate a very simple, two-dimensional optimization problem.  FIGS. 1A–B  show two different views of a two-dimensional surface, with the surface mathematically expressed as follows:
   z=f ( x 1 ,x 2) 
where the domain of the problem is the two-dimensional surface and the goal of the problem is to identify values of x1 and x2 for which the dependent variable z is minimal. In other words, the graph of the surface is a collection of points (x1, x2, z) characterized by a function of two variables, x1 and x2, plotted with respect to the x1-axis  102 , the x2-axis  104 , and the z-axis  106 . As a concrete example, z might represent the cost associated with manufacturing widgets of length x1 and width x2. An optimization problem with regard to this model may be to find the optimal widget length and widget width in order to reduce the cost of manufacturing a widget to the lowest possible value. In most cases, a closed-form expression for the solution, or even a mathematical model, for an optimization problem is not available. Instead, one may begin with some initial, feasible values for the decision variables, in the present case x1 and x2 values corresponding to a point on the two-dimensional surface  110 , and then attempt to find a path, denoted in  FIG. 1A  by the dark curve  112 , from the initial, feasible position  110 , to a local or global minimum value of the function z min  that occurs at a point  114  on the surface with coordinates (x1 min , x2 min , and z min ).
 
   Note that one possible approach for finding a near-optimal or optimal solution to the problem would be to compute z for every possible pair (x1, x2). However, in general, even for extremely simple problems, such an approach would be computationally prohibitive. In many cases, the domains for decision variables, such as decision variables x1 and x2, are subsets of the real numbers, or subsets of other natural domains, which may greatly decrease the potential number of decision-variable-values for which the function would need to be evaluated in order to conduct a complete, combinatory search. However, even in these cases, the combinatory search method is generally not feasible. 
   It should be noted that, although many two-dimensional-surface illustrations, such as  FIGS. 1A–B  are used in the following discussion, they are provided only for illustration purposes. In real-world problems, there may be thousands, hundreds of thousands, or millions of variables and constraints, leading to extremely high-dimensional problem domains. In general, these problem domains may be considered to be hyper-dimensional volumes, or manifolds, bounded by hyper-dimensional surfaces specified by the constraints. It is not possible to illustrate such problem domains, but techniques used to address hyper-dimensional problem domains may analogized to 3-dimensional illustrations. For example, a feasible point in a two-dimensional problem domain is a point of a two-dimensional surface, while a feasible point in a hyper-dimensional problem domain is a point within the hyper-dimensional volume, or manifold, enclosed by the hyper-dimensional boundary surface. 
   Many other approaches are possible for traversing the domain of a mathematical model from an initial, feasible point to a local or global minimum point. One can, for example, compute the gradient for the function at each point, and proceed along the surface in the direction of the gradient in order to most steeply descend along the surface to a local or global minimum. This approach is best suited for convex problems, in which a single global minimum exists. The hyper-dimensional volume domain for a real-word optimization problem is quite often non-convex, and includes many critical points, including saddle points, local maxima, and local minima in addition to a global minimum. Therefore, a steepest-descent approach using gradients may be used to traverse the hyper-dimensional problem domain to a local minimum, but a searching method must nonetheless be employed in order to visit and evaluate multiple local minima. For example, in  FIG. 1A , had the trajectory dictated by the steepest-descent method ended at the local minimum represented by point  116 , then some other approach would be needed in order to detect and evaluate the global minimum at point  112 .  FIG. 1B  shows the problem-domain surface and path to the global minimum shown in  FIG. 1A  rotated by approximately forty-five degrees about the z-axis. 
     FIG. 2  is a control-flow diagram for a highest level view of an optimization method and system that represents one embodiment of the present invention. In a first step,  201 , a particular system or problem to be optimized is mathematically modeled, and the mathematical model serves as input to subsequent optimization steps. Many of the embodiments of the present invention include graphical user interfaces, text-based interpreters, and other such interfaces in order to facilitate mathematical modeling of problems. Below, a particular mathematical model for a hypothetical problem is provided as an example of a mathematical model for an optimization problem. In general, the mathematical model is provided in symbolic form, comprising a set of mathematical expressions including symbolic terms and operators. Next, in step  202 , the method and system undertake a preprocessing step in order to transform the symbolically expressed mathematical model to a standard form needed by methods that represent embodiments of the present invention. Preprocessing steps are discussed below with respect to the hypothetical mathematical model, and the preprocessing steps are formally described. Next, in step  203 , the symbolic standard form of the mathematical model is used to set up an iteration and to iteratively optimize the function obtained using an instance of the interior-point method and a gradient-steered descent-field trajectory in order to identify one or more local minima within a hyper-dimensional problem domain. Various strategies and descent-field mechanisms are employed to find and evaluate a possibly large number of local minima in order to identify a near-optimal solution to the mathematically modeled problem. Various embodiments of the present invention employ different search strategies. It is important to note that neither the methods that represent the embodiments of the present invention, nor any other optimization methods currently employed, are guaranteed to always find a global minimum within a given tolerance or a reasonable amount of time using a reasonable amount of computer resources. However, most practical optimization problems are bounded, in that the number of local minima is finite. Finally, once a best solution is found within some specified tolerance or the maximum amount of time using a specified amount of computing resources, this solution is transformed, in a post-processing step  204 , in order to provide one or more near-optimal solutions in the context of the mathematical model for the problem originally provided in step  201 . 
   The mathematical modeling step and preprocess step,  201  and  202  in  FIG. 2 , respectively, are perhaps best appreciated using a concrete, hypothetical problem.  FIG. 3  shows a hypothetical problem used to illustrate the mathematical modeling and preprocessing steps. The hypothetical problem involves p manufacturing plants,  302 – 304 , w warehouses,  306 – 311 , and c customers,  314 – 318 . A manufacturing plant, such as manufacturing plant  302 , can ship goods to any of the w warehouses, indicated in  FIG. 3  by various types of lines, such as dashed line  320  connecting manufacturing plant  302  and warehouse  306 . Each warehouse, in turn, can ship goods to any of the c customers, with paths between warehouses and customers indicated by various types of lines, such as the dash line  322  interconnecting warehouse  306  and customer  314 . In  FIG. 3 , a different style of line is used to show the path from a particular manufacturing plant to the various warehouses, and from a particular warehouse to the various customers. 
     FIGS. 4–7  illustrate further aspects and characteristics of the hypothetical problem using the illustration conventions of  FIG. 3 . For example, the notation x ij   pw  indicates the amount of goods shipped from manufacturing plant i to warehouse j and terms of the form x jk   wc  indicate the amount of goods shipped from warehouse j to customer k. Although there are w potential warehouses, not all w warehouses need to be used. A Boolean variable b j  indicates whether or not warehouse j is employed, with the Boolean value “1” or “TRUE” indicating that the warehouse is used, and the Boolean value “0” of “FALSE” indicating that the warehouse is not used. In  FIG. 4 , for example, warehouses  1  ( 306 ),  3  ( 308 ), and  5  ( 310 ) are used, as indicated by the solid borders of the squares representing used warehouses and the Boolean values “1” associated with the used warehouses, and warehouses  2  ( 307 ),  4  ( 309 ), and all warehouses following warehouse  5 , including warehouse w  311 , are not used. The used warehouses are alternatively numbered, in  FIG. 4 , with sequential numbers in parentheses so that the three used warehouses are sequentially numbered 1, 2, and n, where n is the total number of used warehouses and equals 3 in the example shown in  FIG. 4 . 
     FIG. 5  illustrates costs associated with shipping goods from manufacturing plants to warehouses and from warehouses to customers. A term of the form c ij   pw  indicates the per item cost of shipping goods from a manufacturing plant i to warehouse j, and a term of the form c jk   wc  indicates the per item cost of shipping goods from warehouse j to customer k. In  FIG. 5 , each of the shipping routes is labeled with the corresponding per-item cost. 
     FIG. 6  illustrates additional characteristics of the hypothetical problem. With each manufacturing plant, there is an associated manufacturing capacity, with a term having the form s i  indicating the shipping capacity of manufacturing plant i. Each warehouse has a receiving capacity, expressed by a term v j  indicating the capacity of warehouse j to receive goods. Also, there is a cost of operation for each warehouse, with a term of the form r j  indicating the cost of operating warehouse j. Finally, there is a demand of goods associated with each customer. A term of the form d k  indicates the demand for goods by customer k. In 
     FIG. 6 , the supply capacities, receiving capacities, cost of operations, and demands associated with manufacturing plants, warehouses, and customers are all explicitly detailed. Finally,  FIG. 7  shows a simple two-manufacturing-plant, three-potential-warehouse, five-customer hypothetical problem associated with numerical costs, capacities, and demands. 
   The optimization problem illustrated in  FIGS. 3–7  is to figure out the lowest cost configuration of operating warehouses, manufacturing-plant-to-warehouse shipping volumes, and warehouse-to-customer shipping volumes such that all the customer demands are met and no capacity is exceeded. A more formal, symbolic mathematical expression of the optimization problem follows. 
   First, the total cost associated with any particular system configuration is expressed as: 
             total   ⁢           ⁢   cost     =         ∑     i   =   1     p     ⁢           ⁢       ∑     j   =   1     n     ⁢           ⁢       c   ij   pw     ⁢     x   ij   pw           +       ∑     j   =   1     n     ⁢           ⁢       ∑     k   =   1     c     ⁢           ⁢       c   jk   wc     ⁢     x   jk   wc           +       ∑     j   =   1     w     ⁢           ⁢       r   j     ⁢     b   j                 
In other words, the total cost is the sum of the costs associated with all shipments from manufacturing plants to warehouses, from warehouses to customers, and the costs of operating warehouses. A minimum is sought in the hypothetical optimization problem for the total cost. In addition, the model includes various constraints. First, a manufacturing plant cannot ship a greater amount of goods than the capacity of the manufacturing plant, expressed as follows:
 
               ∑     j   =   1     n     ⁢           ⁢     x   ij   pw       ≤     s   i           
Customers all receive at least an amount of goods equal to their demands, expressed as follows:
 
               ∑     j   =   1     n     ⁢           ⁢     x   jk   wc       ≥     d   k           
Furthermore, the amount of goods sent by a manufacturing plant to each warehouse is less than or equal to the warehouse&#39;s capacity for receiving goods, expressed as follows:
 
               ∑     i   =   1     p     ⁢     x   ij   pw       ≤     v   j           
Similarly, a warehouse ships an amount of goods less than or equal to the warehouse&#39;s capacity to receive goods, expressed as follows:
 
               ∑     k   =   1     c     ⁢           ⁢     x   jk   wc       ≤     v   j           
Warehouses do not store goods. Each warehouse ships out everything that it receives, expressed by the conservation constraint that follows:
 
                 ∑     i   =   1     p     ⁢           ⁢     x   ij   pw       -       ∑     k   =   1     c     ⁢           ⁢     x   jk   wc         =   0         
The shipment pathways illustrated in  FIGS. 3–7  are unidirectional or, in other words, goods are not returned, expressed as follows:
 x ij   pw ≧0 and  x   jk   wc ≧0 
For the hypothetical problem it is assumed that there is sufficient manufacturing capacity and warehouse-receiving capacity to satisfy the customer demands, expressed as follows:
 
               ∑     i   =   1     p     ⁢           ⁢     s   i       ≥       ∑     j   =   1     w     ⁢           ⁢     v   j       ≥       ∑     k   =   1     c     ⁢           ⁢     d   k             
Goods may be shipped to and from only operating warehouses, expressed by the following two constraints:
   x   ij   pw (1 −b   j )=0   x   jk   wc (1 −b   j )=0 
The number of operating warehouses is less than or equal to the total potential number of warehouses w:
 
               ∑     j   =   1     n     ⁢           ⁢     b   j       ≤   w         
Finally, the terms y j  are Boolean values, expressed as follows:
 b∈{0,1} 
   The above mathematical model may be slightly transformed to more concisely express the cost function and associated constraints. First, the amount-of-goods-shipped parameters of the types x ij   pw  and x jk   wc  may be replaced by scaled shipment terms as follows: 
               m   ij   pw     =       x   ij   pw       s   i         ,     0   ≤     m   ij   pw     ≤   1                     m   jk   wc     =       w   jk   wc       d   k         ,     0   ≤     m   jk   wc     ≤   1           
The warehouse-capacity-based constraints, provided above, turn out to be redundant, and can therefore be removed from the set of constraints without affecting the problem model or solutions obtained for the model. The conservation constraint indicating the warehouses do not store goods is modified to use scaled shipment terms, as follows:
 
               ∑     i   =   1     p     ⁢           ⁢       s   i     ⁢     m   ij   pw         ≥       ∑     k   =   1     c     ⁢           ⁢       d   k     ⁢     m   jk   wc               
Non-operable warehouses do not ship goods, as expressed by the following constraint:
   m   jk   wc   ≦b   j   
The total cost may be expressed, using the scaled shipment of goods terms, as follows:
 
               ∑     i   =   1     p     ⁢           ⁢       ∑     j   =   1     n     ⁢           ⁢       s   i     ⁢     c   ij   pw     ⁢     m   ij   wc           +       ∑     j   =   1     n     ⁢           ⁢       ∑     k   =   1     c     ⁢           ⁢       d   j     ⁢     c   jk   wc     ⁢     m   jk   wc           +       ∑     j   =   1     n     ⁢           ⁢       r   j     ⁢     b   j               
The capacity constraints are expressed, in terms of the scaled shipment of goods terms, as follows:
 
               ∑     j   =   1     n     ⁢     m   ij   pw       ≤   1                   ∑     j   =   1     n     ⁢     m   jk   wc       =   1         
The above-listed capacity constraints can be recast as an inequality with respect to 0, as follows:
 
                 ∑     i   =   1     p     ⁢       s   i     ⁢     r   ij   pw         -       ∑     k   =   1     c     ⁢       d   k     ⁢     m   jk   wc           ≥   0         
The fact that the number of operating warehouses must be less than or equal to the total potential number of warehouses, that goods cannot be returned, and that the b j  terms are Booleans, are now represented as:
 
   
     
       
         
           
             
               ∑ 
               
                 j 
                 = 
                 1 
               
               n 
             
             ⁢ 
             
               b 
               j 
             
           
           ≤ 
           w 
         
       
     
     
       
         
           
             m 
             ij 
             pw 
           
           , 
           
             
               m 
               jk 
               wc 
             
             ≥ 
             0 
           
         
       
     
     
       
         
           
             b 
             j 
           
           ∈ 
           
             { 
             
               0 
               , 
               1 
             
             } 
           
         
       
     
   
     FIG. 8  is a more detailed control-flow diagram for the preprocess routine invoked in step  202  of the high-level control-flow diagram shown in  FIG. 2 . Each step of the preprocess routine illustrated in  FIG. 8  is further discussed below with reference to corresponding transformations of the hypothetical mathematical model described above. First, in step  801 , the preprocess routine receives the initial, symbolic, mathematically expressed model for the optimization problem, including constraints. In step  803 , the preprocessing routine attempts to remove as many redundant constraints as possible. An example of the removal of redundant constraints occurs in the above transformed model, with the warehouse-capacity constraints in the original form of the model recognized as redundant and removed from the transformed model. Next, in step  805 , all Boolean and integer-value decision variables, such as the Boolean variables b j , in the hypothetical example, are transformed into floating-point variables. In step  807 , the resulting symbolic, mathematically expressed model is placed in a standard form. In step  809 , equality constraints are replaced by relaxed inequality constraints. Each equality constraint is relaxed with three inequality constraints using an additional positively valued variable. In step  811 , an additional variable z is added, along with an additional constraint to convert the criterion into an additional constraint. Finally, in step  813 , the transformed optimization problem is converted into an Interior Point Method barrier function by adding together terms. Each term is an evaluation of a barrier function in terms of an expression representing a constraint multiplied by a control variable. This Interior Point Method Barrier function contains two sets of control variables u and r. As a result of the preprocessing routine, the mathematical model is symbolically expressed in a form that can be optimized by the iterative optimization technique shown as step  203  in  FIG. 2 . 
     FIG. 9  shows the initial symbolic, mathematically expressed model for the specific optimization problem illustrated in  FIG. 7 . The cost function,  902 , is to be optimized with respect to vectors of integer variables x and Boolean variables b. The constraints include manufacturing-plant-capacity constraints  904 , customer-demand constraints  906 , conservation constraints  908 , constraints  910  that model the fact that no flow can occur to or from a warehouse that is not selected by the optimization, the corresponding Boolean variable being zero, a constraint that only two of the three warehouses may be operational  912 , a constraint specifying that there are no return shipments of goods from customers to warehouses or from warehouses to manufacturing plants  914 , and the constraint that the b j  are Boolean-valued variables  916 . 
   The third step of the preprocessing routine, step  805  in  FIG. 8 , involves transforming Boolean, integer, and other non-floating-point variables to floating-point variables. Each Boolean variable b i  is replaced with a floating-point value y i  in the range [0,1] and with an additional constraint:
 
 y   i (1 −y   i )=0
 
The additional constraint implies that the floating-point variable y i  can only have one of the two values 0 or 1. Each integer variable is replaced by an equivalent base-2 polynomial with k terms, where k is approximately equal to the base-2 logarithm of the maximum value that the integer variable can have. Each term of the base-2 polynomial comprises a power of 2 multiplied by a Boolean coefficient, and thus the base-2 polynomial is, in turn, replaced with k floating-point variables y 0 , y 1 , . . . y k−1 , all within the range [0,1] along with k constraints in the form:
 
 y   i (1 −y   i )=0
 
The conversion of decision variables of various types to floating-point numbers is summarized formally as follows:
 
             x   i     ∈     R   ⇒       x   i     ⁢           ⁢   is   ⁢           ⁢   float                     x   i     ∉     R   →       x   i     ⁢           ⁢   is   ⁢           ⁢   transformed   ⁢           ⁢   to   ⁢           ⁢   float   ⁢     :                       ∀       x   i     ∈   Bool       ,         x   i     ⇒     y   i       ∈     [     0   ,   1     ]       ,         y   i     ⁡     (     1   -     y   i       )       =   0                   ∀       x   i     ∈     Z   +         ,         x   i     ⇒       ∑     i   =   0     k     ⁢       b   i     ⁢     2   i           =         b   0     ⁢     2   0       +       b   1     ⁢     2   1       +       b   2     ⁢     2   2       +   …   +       b     k   -   1       ⁢     2     k   -   1                           where   ⁢           ⁢     maxvalue   ⁡     (     x   i     )         =       2   k     -   1                         ⁢         log   2     ⁡     (     maxvalue   ⁡     (     x   i     )       )       ≅   k                         ⁢       y   0     ,     y   1     ,           y   2     ⁢           ⁢   …   ⁢           ⁢     y     k   -   1         ∈       [     0   ,   1     ]     ⁢           ⁢   and   ⁢           ⁢     x   i         =       ∑     j   =   0       k   -   1       ⁢       y   j     ⁢     2   j                                 ⁢         y   0     ⁡     (     1   -     y   0       )       =   0                         ⁢         y   1     ⁡     (     1   -     y   1       )       =   0                         ⁢   ⋮                       ⁢         y     k   -   1       ⁡     (     1   -     y     k   -   1         )       =   0                   x   i     ∈     R   ⇒       x   i     ⁢           ⁢   is   ⁢           ⁢   float                     x   i     ∉     R   →       x   i     ⁢           ⁢   is   ⁢           ⁢   transformed   ⁢           ⁢   to   ⁢           ⁢   floats   ⁢     :                       ∀       x   i     ∈   Bool       ,         x   i     ⇒     y   i       ∈     [     0   ,   1     ]       ,         y   i     ⁡     (     1   -     y   i       )       =   0                   ∀       x   i     ∈     Z   +         ,         x   i     ⇒       ∑     i   =   0     k     ⁢       b   i     ⁢     2   i           =         b   0     ⁢     2   0       +       b   1     ⁢     2   1       +       b   2     ⁢     2   2       +   …   +       b     k   -   1       ⁢     2     k   -   1                           where   ⁢           ⁢     maxvalue   ⁡     (     x   i     )         =       2   k     -   1                         ⁢         log   2     ⁡     (     maxvalue   ⁡     (     x   i     )       )       ≅   k                         ⁢       y   0     ,     y   1     ,           y   2     ⁢   …   ⁢           ⁢     y     k   -   1         ∈       [     0   ,   1     ]     ⁢           ⁢   and   ⁢           ⁢     x   i         =       ∑     j   =   0       k   -   1       ⁢       y   j     ⁢     2   j                                 ⁢         y   0     ⁡     (     1   -     y   0       )       =   0                         ⁢         y   1     ⁡     (     1   -     y   1       )       =   0                         ⁢   ⋮                       ⁢         y     k   -   1       ⁡     (     1   -     y     k   -   1         )       =   0             FIG. 10  shows transformation of the hypothetical mathematical model into a model employing only floating-point variables.
 
   Next, the mathematical model is placed in standard form. The formal definition for a standard-form symbolic problem representation is provided as follows:
         minimize C(x)   subject to:       

                   ⁢             g   ~     i     ⁡     (   x   )       ≥     0   ⁢           ⁢   for   ⁢           ⁢   i       =       0   ⁢           ⁢   to   ⁢           ⁢     M             ⁢   1         -   1                                 ⁢     h   ~       i     ⁢     (   x   )       =       0   ⁢           ⁢   for   ⁢           ⁢   i     =         M             ⁢   1       ⁢           ⁢   to   ⁢           ⁢     M             ⁢   2         -   1                 where             x   =         [                   ⁢     x             ⁢   0                           ⁢     x             ⁢   0                           ⁢     x             ⁢   0                 ⋮                     ⁢     x             ⁢   0               ]     ⁢           ⁢   and   ⁢           ⁢     x             ⁢   0         ∈   ℛ           
In other words, a standard form model includes a cost function C(x), a number M 1  of inequality constraints of the form
   {tilde over (g)}   i ( x )≧0 
and a number of equality constraints of the form:
   {tilde over (h)}   i ( x )=0   FIG. 11  shows the mathematical model for the hypothetical problem transformed to standard form. Note that transformation of the original equality and inequality constraints to constraints with 0 on the right-hand side, as needed for the standard form, is a simple matter of adding terms to both sides of the constraint equations that are not already in standard form.
 
   The next step in the preprocessing routine ( 809  in  FIG. 8 ) is to change all equality constraints to inequality constraints. In order to do so, each equality constraint is transformed to two inequality constraints, each inequality constraint involving the left-hand expression of the original equality constraint and a positive binding variable. Transformation of an equality constraint h j (x) to two inequality constraints is formally expressed as follows: 
               h   j     ⁡     (   x   )       →       [                 h   j     ⁡     (   x   )       +     s   j       ≥   0                   -       h   j     ⁡     (   x   )         +     s   j       ≥   0           ]     →     [             g     m   +   j       ⁡     (     x   ,     s   j       )                   g     m   +   j   +   1       ⁡     (     x   ,   ,     s   j       )             ]             
In addition for each j, the constraint
 s j ≧0 
is added to ensure that the constraints expressing h j (x) are feasible.
 
   Thus, following transformation of equality constraints to inequality constraints, the total number of constraints increases by a number 2k equal to the original number of equality constraints in the standard-form model of the problem.  FIGS. 12A–D  illustrate the effect of binding variables on transformed equality constraints.  FIG. 12A  shows the original inequality constraint h j (x) with the value 0 ( 1202 ). When the equality constraint h j (x) is transformed to three inequality constraints, as discussed above, then the range of values for h j (x) expands about 0 to include all values between −s j    1204  and s j    1206 . Thus, transforming an equality into three inequality constraints essentially relaxes the equality constraint from a point to an interval. Relaxing equality constraints facilitates the iterative optimization methods, to be described below, by preventing equality constraints from overconstraining, or dominating, the initial descent towards one or more local minima. Relaxing equality constraints also facilitates an interior-point descent, since equality constraints force optimization trajectories to problem-domain boundaries. However, as the optimization begins to converge, the relaxation of equality constraints needs to be removed. This is accomplished, as shown in  FIGS. 12C and 12D , by decreasing the absolute value of the barrier variable sj in order to more and more tightly constrain the value of h j (x) to a small region symmetrically disposed about 0. The degree of relaxation of an equality constraint can be expressed by a coefficient r i  for each binding variable s j  added to the criterion. Transformation of equality constraints to inequality constraints can be formally expressed as:
         minimize C(x)   constrained by:       

                   g   i     ⁡     (   x   )       ≥     0   ⁢           ⁢   i       =   0     ,   1   ,   …   ⁢           ,     m   -   1                       h   j     ⁡     (   x   )       =       0   ⁢           ⁢   j     =   m       ,   1   ,   …   ⁢           ,       m   +   k   -       1   ⁢     
     ⁢           ⇓     
     ⁢   minimize     ⁢           ⁢       C   ~     ⁡     (   x   )           =       C   ⁡     (   x   )       +       ∑     i   =   0         2   ⁢   k     -   1       ⁢       r   i     ⁢     s   i                           constrained   ⁢           ⁢   by   ⁢     :       ⁢                               ⁢             g   i     ⁡     (   x   )       ≥     0   ⁢           ⁢   i       =   0     ,   1   ,   …   ⁢           ,     m   +     2   ⁢   k     -   1               FIG. 13  shows the symbolic, mathematically expressed model for the hypothetical problem following transformation of equality constraints to inequality constraints.  FIG. 14  graphically illustrates the effect of the binding variables r i . In  FIG. 14 , the optimization method begins at initial point  1402  and descends through a problem-domain  1404  bounded by the cup-shaped surface  1406  towards a minimum  1408 . Initially, the coefficients r i  for the binding variables s i  are small, providing a large total range of values, shown in  FIG. 14  by the disk-shaped region  1410 , about the initial point  1402 . As the descent through the problem domain approaches the minimum  1408 , the coefficients r i  of the binding variables se are increased, more and more tightly constraining the solution sets until, near the optimum  1408 , the transformed equality constraints again become extremely tightly constrained  1412  because as the optimization method proceeds the s i  become smaller and smaller.
 
   The next step of the preprocessing routine involves adding an envelope variable z along with an additional constraint to the mathematical model of the problem ( 811  in  FIG. 8 ). The envelope variable z is defined as follows:
 
 z≧{tilde over (C)} ( x )
 
This leads to the following inequality constraint:
 
 z−{tilde over (C)} ( x )≧0
 
Following introduction of the envelope variable z, the formal statement of the problem becomes:
         minimize z   with respect to x 0 , x 1  . . . , x n−1 , s 0 , s 1 , . . . , s k−1 , z   constrained by:
 
 g   i ( x ) i=0 to m−1
 
 g   j ( x,s   j ) j=0 to 2k−1
 
 z−{tilde over (c)} ( x )≧0
 
 FIG. 15  shows the mathematical model for the hypothetical problem following introduction of the envelope variable z.
       

     FIGS. 16A–D  illustrate the concept of barrier functions. Assume, for example, that the problem domain for an optimization problem is the wedge-shaped volume with a parabolic base  1602  shown on  FIG. 16A . An initial point for an optimization  1604  has a relatively large z-value  1601 , which decreases to 0 at the minimum  1608 . In  FIG. 16B  shows an x-y plane projection  1610  of the problem domain  1602 . The optimization trajectory from the initial point  1604  to the minimal point  1608  is shown as three line segments  1612 – 1614 . In an interior point method, when traversing the problem domain, the optimization trajectory must remain within the problem domain. In other words, the current point in a descent trajectory through the problem domain to a minimum must always remain within the problem domain, interior to the boundary of the manifold in a hyper-dimensional problem domain. For example, as shown in  FIG. 16B , the initial trajectory  1612  is angled towards the edge of the volume corresponding to the problem domain, and would intersect the boundary at point  1616  and, if continued in the same direction, would pass through the boundary to points exterior to the problem domain along the dashed line segment  1618 . In order to keep the trajectory of the optimization within the problem domain, it is necessary to include barrier functions within a problem model that have little or no effect within the problem domain, but rise to large values as the optimization trajectory nears the boundary surface of the problem domain, in order to deflect the trajectory back from the boundary into the problem domain. For example, in  FIG. 16B , as the optimization trajectory  1612  approaches point  1616 , the barrier functions added to the mathematical model begin to rapidly increase in order to drive the value of the envelope variable z to high values and therefore deflect the optimization trajectory away from the inner surface of the problem domain. Following deflection, a new optimization trajectory  1613  is established which leads, through continued optimization, to a new second approach to a boundary point  1620 , just prior to which the barrier constraints again force the value of the objective function z to large values in order to again deflect the optimization trajectory back into the interior of the problem domain and to a new, final trajectory  1614 . 
     FIGS. 16C–D  illustrate the operation of the barrier functions. Assuming a current optimization trajectory  1622 , as shown in  FIG. 16C , consider the value of the objective function z with respect to the x-coordinate, shown in  FIG. 16D , for a traversal of the problem of the domain along the current trajectory  1622  shown in  FIG. 16C . As can be seen in  FIG. 16D , the value of z is relatively low for values of x within the interior of the problem domain along the current trajectory  1622 . However, as x decreases towards the left boundary point  1624  or increases towards the right boundary point  1626 , the barrier functions have the effect of rapidly increasing the value of z  1628  and  1639 , respectively. In other words, for each inequality in the model, a term involving a barrier function is added to the interior point objective function in order to constrain trajectories to the interior of the problem domain. If a trajectory leaves the problem domain, or reaches the exterior of the enclosing surface of the problem domain, the optimization methods that represent some of the embodiments of the present invention could fail. 
   Adding terms involving barrier functions to the model is a simple way to ensure that the optimization-trajectory stays within the problem domain. A barrier function is added for each constraint, since it is the constraints that specify the surface enclosing the problem domain. A suitable barrier-function term for the optimization techniques employed in various embodiments of the present invention is the negative of a positive variable coefficient u i , the control variable associated with the term, times the logarithm of the constraint expression g i , or, in other words:
 
−u i  log(g i )
 
As the value of the constraint expression g i  decreases towards 0, the value of the logarithm of the constraint log(g i ) decreases rapidly with a corresponding rapid increase in absolute value. By taking the negative of the log (g i ), a rapidly increasing value is obtained as the constraint expression g i  closely approaches 0, approaching infinity in the limit, thus preventing the constraint from ever decreasing to 0. The variable coefficients u i  can be considered to be control variables, much like the coefficients r i  of the above-described binding variables. The u i  and r i  coefficients are modified, during an optimization, in order to steer the optimization trajectory in a most efficient manner towards local minima. However, unlike the r i  coefficients, the u i  coefficients are initially large, and are decreased by one of the embodiments of this invention. as the optimization trajectory approaches a local minimum, since local minima often fall close to, or on the surface of, the problem domain. A formal expression the addition of barrier functions to the objective function is provided as follows:
 
               minimize     x   ,   s   ,   z       ⁢     F   ⁡     (     x   ,   s   ,   z   ,   u   ,   r     )         =     z   -       ∑     i   =   0       i   =     m   +     2   ⁢   k           ⁢       u   i     ⁢           ⁢     log   ⁡     (       g   i     ⁡     (     x   ,   s   ,   z   ,   r     )       )                     FIG. 17  shows the symbolic, mathematical model for the hypothetical problem following addition of barrier constraints. The addition of the barrier functions completes the preprocessing steps embodied in the preprocessing routine in a first, described embodiment.
 
   Next, an overview of the core optimization method that represents one embodiment of the present invention is provided.  FIG. 18  illustrates the use of a global gradient descent field within a problem domain in order to provide a direction for an optimization trajectory leading from an initial point to a local minimum. At each point within the problem domain, such as point  1802 , the partial derivative vector of the objective function F(x,s,z,u,r) may be calculated as a vector of first partial derivatives of the objective function with respect to the decision variables x i , the binding variables s i , and the envelope variable z. Similarly the Hessian of the function F(x,s,z,u,r) given as the matrix of second partial derivatives of F(x,s,z,u,r) with respect to the decision variables x i , the binding variables s i , and the envelope variable z can be computed at every point of the domain The gradient of the objective function is defined as the vector at each point in the domain computed as the product of the inverse of the Hessian and the partial derivative of F(x,s,z,u,r). It provides a direction  1804  from the point  1802  representing a descent direction in the value of the objective function with respect to the decision variables and binding variables. Since a gradient can be calculated for each point within the problem domain, the gradient of the objective function specifies a vector field, pointing in a direction of descent and, specifying a path to a minimum within the problem domain. In general, gradients for points, such as point  1802  in  FIG. 18 , far from local minima have relatively large absolute values, while gradients for points close to a local minima, such as point  1806 , have gradients with relatively small absolute values. 
   The mathematical description of the optimization methods is facilitated by a number of notational conventions. First, the objective function
         F(x,s,z,u,r)
 
can be expressed as:
       

           f   ⁡     (     y   ,   u   ,   r     )                   with   ⁢     
     [         x           s           z         ]     =       y   ⁢     
     [           u   0               u   1             ⋮             u     m   +     2   ⁢   k               ]     =       u   ⁢     
     [           r   0               r   1             ⋮             r       2   ⁢   k     -   1             ]     =   r             
Thus, the objective function can be thought of as a function of a vector y of decision variables, binding variables, and the envelope variable, the vector u, a vector of barrier-function coefficients, and the vector r, a vector of binding variable coefficients. The vector y is referred to as the “state vector,” and the vectors u and r are referred to as “control vectors.” The vector y has N elements, where N=n (decision variables)+k (binding variables)+1 (envelope variable).
 
     FIG. 19  illustrates an optimization trajectory. The trajectory  1902  is plotted with respect to the lengths of vectors r and u as horizontal axes, and the value of the objective function F as a vertical axis. The optimization trajectory begins at an initial point  1904  and descends in the value of the objective function F towards a local minimum. During the descent, the value of the barrier functions expressed by the length of vector u, decreases towards 0, while the value of the coefficients of the binding variables, as expressed as the length of vector r, increases. The trajectory is divided into segments of length Δτ, such as a first segment  1906 . During optimization, the vectors y and u are continuously recalculated to drive the value of the objective function lower. In  FIG. 19 , the value of the binding variable coefficients, represented as vector r, is also shown to be recalculated during each Δτ interval, although in the described embodiment of the present invention, the binding variable coefficients are adjusted only towards the end of an optimization trajectory. 
   The discrete optimization technique may be expressed as follows:
 
 y   k+1   =y   k +φ k ( y   k   , r   k   , u   k )
 
where k=0, 1, . . . , k convergence . In other words, a new state-vector value y k+1  is computed, in each iteration of a computational optimization method, from a current state-vector value y k  and a discrete function φ k  that depends on y k , r k , and u k , with the iteration count k ranging from 0 to a value k convergence  at which point the state vector for a local minimum of the objective function might be determined. In a first step of optimization, the state and control vectors, y, u, and r may be parameterized by a continuous iteration parameter τ, where τ ranges from 0 to an iteration horizon T. By parameterizing vectors y, U, and r by the continuous iteration variable τ, the iteration above is transformed to a first order differential equation in which the y variable is the dependent variable and u, and r are control variables. Well-known differential-equation-based mathematical technique can be used for computing efficient optimization trajectories as solutions of this differential equation. This is discussed below. The objective function is now expressed formally as:
 
 F ( y (τ),  u (τ),  r (τ))
 
where 0≦τ≦T. For example if iteration y k+1 =y k   +φ   k (y k , r k , u k ) is the Newton iteration, the corresponding differential equation has the following form
 
   
     
       
         
           
             
               ⅆ 
               y 
             
             
               ⅆ 
               τ 
             
           
           = 
           
             
               
                 H 
                 
                   - 
                   1 
                 
               
               ⁡ 
               
                 ( 
                 
                   
                     y 
                     ⁡ 
                     
                       ( 
                       τ 
                       ) 
                     
                   
                   , 
                   
                     u 
                     ⁡ 
                     
                       ( 
                       τ 
                       ) 
                     
                   
                   , 
                   
                     r 
                     ⁡ 
                     
                       ( 
                       τ 
                       ) 
                     
                   
                 
                 ) 
               
             
             ⁢ 
             
               G 
               ⁡ 
               
                 ( 
                 
                   
                     y 
                     ⁡ 
                     
                       ( 
                       τ 
                       ) 
                     
                   
                   , 
                   
                     u 
                     ⁡ 
                     
                       ( 
                       τ 
                       ) 
                     
                   
                   , 
                   
                     r 
                     ⁡ 
                     
                       ( 
                       τ 
                       ) 
                     
                   
                 
                 ) 
               
             
           
         
       
     
   
   It should be noted that Newton&#39;s method is but one approach to construct the differential equation providing an expression for 
               ⅆ   y       ⅆ   τ       .         
It was chosen for simplicity and for ease of computation. However, any of many other approaches may be used in embodiments of the present invention for approximating
 
   
     
       
         
           
             
               ⅆ 
               y 
             
             
               ⅆ 
               τ 
             
           
           . 
         
       
     
   
   In the following discussion, the partial derivative vector G and Hessian matrix H for the objective function F are frequently used. The partial derivative vector G is defined as follows: 
                   G   ⁢           ⁢     (       y   ⁢           ⁢     (   τ   )       ,     u   ⁢           ⁢     (   τ   )       ,     r   ⁢           ⁢     (   τ   )         )       =       (         ∂   F     ⁢           ⁢     (       y   ⁢           ⁢     (   τ   )       ,     u   ⁢           ⁢     (   τ   )       ,     r   ⁢           ⁢     (   τ   )         )           ∂   y     ⁢           ⁢     (   τ   )         )     T                 =     [               ∂   F     ⁢           ⁢     (       y   ⁢           ⁢     (   τ   )       ,     u   ⁢           ⁢     (   τ   )       ,     r   ⁢           ⁢     (   τ   )         )         ∂       y   0     ⁡     (   τ   )                         ∂   F     ⁢           ⁢     (       y   ⁢           ⁢     (   τ   )       ,     u   ⁢           ⁢     (   τ   )       ,     r   ⁢           ⁢     (   τ   )         )         ∂       y   1     ⁡     (   τ   )                   ⋮                 ∂   F     ⁢           ⁢     (       y   ⁢           ⁢     (   τ   )       ,     u   ⁢           ⁢     (   τ   )       ,     r   ⁢           ⁢     (   τ   )         )         ∂       y     N   -   1       ⁡     (   τ   )                 ]                 
The Hessian matrix H is a matrix of size N×N, defined as follows:
 
                   H   ⁢     (       y   ⁢     (   τ   )       ,     u   ⁢     (   τ   )       ,     r   ⁢     (   τ   )         )       =         [           ⁢       ∂   G       ∂   y       ]     T     ⁢     (       y   ⁢     (   τ   )       ,     u   ⁢     (   τ   )       ,     r   ⁢     (   τ   )         )                     =     [                   ⁢       ∂     G   0         ∂     y   0                         ⁢       ∂     G   1         ∂     y   0                         ⁢       ∂     G   2         ∂     y   0               ⋯                 ⁢       ∂     G     N   ⁢           -           ⁢   1           ∂     y   0                             ⁢       ∂     G   0         ∂     y   1                         ⁢       ∂     G   1         ∂     y   1                         ⁢       ∂     G   2         ∂     y   1               ⋯                 ⁢       ∂     G     N   ⁢           -           ⁢   1           ∂     y   1                   ⋮       ⋮       ⋮                   ⋮                     ⁢       ∂     G   0         ∂     y     N   -   1                           ⁢       ∂     G   1         ∂     y     N   -   1                           ⁢       ∂     G   2         ∂     y     N   -   1                 ⋯                 ⁢       ∂     G     N   ⁢           -           ⁢   1           ∂     y     N   -   1                   ]       ⁢           ⁢     
     ⁢           ⁢     (       y   ⁢     (   τ   )       ,     u   ⁢     (   τ   )       ,     r   ⁢     (   τ   )         )                   =       [             ∂     G   0         ∂     y   0                 ∂     G   0         ∂     y   1             ⋯       ⋯           ∂     G   0         ∂     y     N   -   1                       ∂     G   1         ∂     y   0                 ∂     G   1         ∂     y   1             ⋯       ⋯           ∂     G   1         ∂     y     N   -   1                       ∂     G   2         ∂     y   0                 ∂     G   2         ∂     y   1             ⋯       ⋯           ∂     G   2         ∂     y     N   -   1                   ⋮       ⋮                               ⋮               ∂     G     N   ⁢           -           ⁢   1           ∂     y   0                 ∂     G     N   -   1           ∂     y   1             ⋯       ⋯           ∂     G     N   -   1           ∂     y     N   -   1                 ]     T       ⁢     
     ⁢     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )       ,     r   ⁡     (   τ   )         )                 =     [               ∂   2     ⁢     F   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )       ,     r   ⁡     (   τ   )         )           ∂     y   0   2             ⋯             ∂   2     ⁢     F   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )       ,     r   ⁡     (   τ   )         )             ∂     y   0       ⁢     ∂     y     N   -   1                     ⋮                   ⋮           ⋮                   ⋮                 ∂   2     ⁢     F   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )       ,     r   ⁡     (   τ   )         )             ∂     y     N   -   1         ⁢     ∂     y   0               ⋯             ∂   2     ⁢     F   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )       ,     r   ⁡     (   τ   )         )             ∂     y     N   -   1         ⁢     ∂     y     N   -   1                   ]                 
The normalized rate of change of the state vector y with respect to the continuous iteration variable τ,
 
               ⅆ   y       ⅆ   τ       =       -         H     -   1       ⁢   G         G   T     ⁢     H     -   1       ⁢   G         ⁢     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )       ,     r   ⁡     (   τ   )         )             
when expressed in terms of the Hessian matrix H and vector G, or as:
 
   
     
       
         
           
             
               ⅆ 
               y 
             
             
               ⅆ 
               τ 
             
           
           = 
           
             - 
             
               
                 w 
                 ⁡ 
                 
                   ( 
                   τ 
                   ) 
                 
               
               
                 ρ 
                 ⁡ 
                 
                   ( 
                   
                     
                       y 
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       u 
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       r 
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                   
                   ) 
                 
               
             
           
         
       
     
     
       
         
           
             where 
             ⁢ 
             
                 
             
             ⁢ 
             
               w 
               ⁡ 
               
                 ( 
                 τ 
                 ) 
               
             
           
           = 
           
             
               H 
               
                 - 
                 1 
               
             
             ⁢ 
             
               G 
               ⁡ 
               
                 ( 
                 
                   
                     y 
                     ⁡ 
                     
                       ( 
                       τ 
                       ) 
                     
                   
                   , 
                   
                     u 
                     ⁡ 
                     
                       ( 
                       τ 
                       ) 
                     
                   
                   , 
                   
                     r 
                     ⁡ 
                     
                       ( 
                       τ 
                       ) 
                     
                   
                 
                 ) 
               
             
             ⁢ 
             
                 
             
             ⁢ 
             and 
           
         
       
     
     
       
         
           
             ρ 
             ⁡ 
             
               ( 
               
                 
                   y 
                   ⁡ 
                   
                     ( 
                     τ 
                     ) 
                   
                 
                 , 
                 
                   u 
                   ⁡ 
                   
                     ( 
                     τ 
                     ) 
                   
                 
                 , 
                 
                   r 
                   ⁡ 
                   
                     ( 
                     τ 
                     ) 
                   
                 
               
               ) 
             
           
           = 
           
             
               G 
               T 
             
             ⁢ 
             
               H 
               
                 - 
                 1 
               
             
             ⁢ 
             
               
                 G 
                 ⁡ 
                 
                   ( 
                   
                     
                       y 
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       u 
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       r 
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                   
                   ) 
                 
               
               . 
             
           
         
       
     
   
   The optimization problem can be recast into a minimum time problem which is a special case of a generic, optimal control problem. In other words, by continualization of the objective function F(y,U,r), continuous methods can be used to obtain a computational approach for computing an optimization trajectory y(τ) consisting of state-variable changes that bring the value of the objective function F(y,u,r) from an initial point to a minimum. However, the problem domain is enormous, and hyper-dimensional, and so, without further computational constraints, or meta-level constraints, although a method for descending the gradient field within the problem domain has been obtained, the problem of computationally descending the gradient field is slow convergent. 
   An optimal control problem is formulated. It seeks to speed-up the convergence by minimizing the convergence horizon and satisfying the requirements on the control variables. In order to control minimization of the objective function F(y(τ), u(τ), r(τ)) the optimal control formulation is used to steer, using the control variables u(τ), r(τ), the optimization trajectory by selecting values for the control variables that minimize the convergence horizon. The optimal control problem for controlling the optimization trajectory is formally expressed as follows: 
               min     u   ,   r       ⁢     ??   ⁡     (       u   ⁡     (   τ   )       ,     r   ⁡     (   τ   )       ,   T     )         =         ∫   o   T     ⁢     ⅆ   τ       +       1   2     ⁢     G   T     ⁢     QG   ⁡     (       y   ⁡     (   T   )       ,     u   ⁡     (   T   )       ,     r   ⁡     (   T   )         )                 
with respect to u(τ), r(τ), τ∈[0, T], where, Q is a symmetric positive definite matrix, subject to the constraints:
 
                     ∂   y       ∂   τ       =         -       w   ⁡     (   τ   )         ρ   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )       ,     r   ⁡     (   τ   )         )           ⁢           ⁢   for   ⁢           ⁢   τ     ∈     [     0   ,   T     ]               (   1   )               u i ∈(0, u max     i   ) ∀i  (2)   u   i (τ+Δτ)≦ u   i (τ) ∀ i   (3) r i ∈[r min     i   , ∞) ∀i  (4)   r   i (τ+Δτ)≧ r   i (τ) ∀ i   (5) 
Minimization of the objective functional for the minimum convergence horizon problem,  (u(τ), r(τ), T), forces F(y(τ), u(τ), r(τ)) to converge to a minimum, where G=0, in the smallest computational horizon T, where the convergence interval is the integral of the continuous iteration variable. The solution to the optimal control problem provides an optimal control policy u*(τ) and r*(τ), τ∈[0, T] that generates a state-vector-change optimization trajectory y*(τ). The first term in the integral defining the objective function  (u(τ), r(τ), T) forces the optimization trajectory to an efficient as possible trajectory y* with respect to the continuous iteration variable τ. The second term within the integral forces the vector G towards 0, thereby forcing the objective function F to a minimum. The combination of these terms therefore produces an optimal control policy u*(τ) and r*(τ) to force the objective function F to a local minimum through an optimal state-vector-change trajectory y*(τ). Now the constraints of the optimal control problem (1)–(5) are discussed. The first constraint involves the normalized steepest descent for the state-vector-change trajectory, and is the primary vehicle by which the original optimization problem is incorporated into the new minimum-time control problem. The second and third constraints ensure that the control vector values u remains non-negative and less than some maximal control value u max     i    and that the u-vector values do not increase during the optimization process. The fourth and fifth constraints are analogous to the second interior constraints with respect to the binding variable coefficients r i . The fourth and fifth constraints ensure that the binding variable coefficients r i  do not fall below a minimum value r min     i   , and that the binding variable coefficients r i  do not decrease during the optimization process.
 
   A primary result of optimal control theory, known as the Pontryagin “Minimum Principle,” provides a framework for the minimum convergence horizon problem to produce an optimal control policy u*(τ) and r*(τ)τ∈[0, T] for finding an optimal state-vector-change optimization trajectory y*(r). This approach involves the Hamiltonian function, defined as follows: 
           ℋ   =     1   -         p   T     ⁡     (   τ   )       ⁢       w   ⁡     (   τ   )         ρ   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )       ,     r   ⁡     (   τ   )         )                   
The Hamiltonian has a minimum value as a function of u and r at u*(τ) and r*(τ). In other words, the Hamiltonian has a minimum value for the optimal control policy. When the control policy u*(τ) is optimal, there exists an absolutely continuous function p(τ), called the costate, such that the derivatives of the state and costate functions y(τ) and p(τ) are as follows for all τ∈[0, T]:
 
               ⅆ       y   *     ⁡     (   τ   )           ⅆ   τ       =       ∂     ℋ   ⁡     (         y   *     ⁡     (   τ   )       ,       u   *     ⁡     (   τ   )       ,     p   ⁡     (   τ   )         )           ∂   p                       ⅆ     p   ⁡     (   τ   )           ⅆ   τ       =     -     (       ∂     ℋ   ⁡     (         y   *     ⁡     (   τ   )       ,       u   *     ⁡     (   τ   )       ,     p   ⁡     (   τ   )         )           ∂   y       )             
satisfying the boundary conditions:
 
               y   *     ⁡     (   0   )       =     y             ⁢   0                     p   ⁢     (   T   )       =         ∂           ⁢     G             ⁢   T         ⁢     (       y   ⁢     (   T   )       ,           ⁢     u   ⁢     (   T   )       ,           ⁢     r   ⁢     (   T   )         )     ⁢           ⁢   QG   ⁢     (       y   ⁢     (   T   )       ,           ⁢     u   ⁢     (   T   )       ,           ⁢     r   ⁢     (   T   )         )         ∂   y                     for   ⁢           ⁢   all   ⁢             ⁢             ⁢   u   ⁢     (   τ   )       ,       ℋ   ⁢     (         y   *     ⁢     (   τ   )       ,       u   *     ⁢     (   τ   )       ,     p   ⁢     (   τ   )         )       ≤       ℋ   ⁡     (         y   *     ⁢     (   τ   )       ,     u   ⁢     (   τ   )       ,     p   ⁢     (   τ   )         )       .             
The control problem solution is summarized as follows:
 
   
     
       
         
           
             
               ⅆ 
               
                 y 
                 * 
               
             
             
               ⅆ 
               τ 
             
           
           = 
           
             - 
             
               w 
               ρ 
             
           
         
       
     
     
       
         
           
             
               ⅆ 
               p 
             
             
               ⅆ 
               τ 
             
           
           = 
           
             
               ∂ 
               
                 ∂ 
                 y 
               
             
             ⁢ 
             
               ( 
               
                 w 
                 ρ 
               
               ) 
             
             ⁢ 
             
               p 
               ⁡ 
               
                 ( 
                 τ 
                 ) 
               
             
           
         
       
     
     
       
         
           
             
               y 
               * 
             
             ⁡ 
             
               ( 
               0 
               ) 
             
           
           = 
           
             y 
             0 
           
         
       
     
     
       
         
           
             p 
             ⁡ 
             
               ( 
               T 
               ) 
             
           
           = 
           
             
               H 
               ⁡ 
               
                 ( 
                 T 
                 ) 
               
             
             ⁢ 
             
               QG 
               ⁡ 
               
                 ( 
                 T 
                 ) 
               
             
           
         
       
     
     
       
         
           
             
               p 
               T 
             
             ⁢ 
             
               
                 w 
                 ⁡ 
                 
                   ( 
                   
                     
                       
                         y 
                         * 
                       
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       
                         u 
                         * 
                       
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       
                         r 
                         * 
                       
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                   
                   ) 
                 
               
               
                 ρ 
                 ⁡ 
                 
                   ( 
                   
                     
                       
                         y 
                         * 
                       
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       
                         u 
                         * 
                       
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       
                         r 
                         * 
                       
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                   
                   ) 
                 
               
             
           
           ≥ 
           
             
               p 
               T 
             
             ⁢ 
             
               
                 w 
                 ⁡ 
                 
                   ( 
                   
                     
                       
                         y 
                         * 
                       
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       u 
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       r 
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                   
                   ) 
                 
               
               
                 ρ 
                 ⁡ 
                 
                   ( 
                   
                     
                       
                         y 
                         * 
                       
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       u 
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                     , 
                     
                       r 
                       ⁡ 
                       
                         ( 
                         τ 
                         ) 
                       
                     
                   
                   ) 
                 
               
             
           
         
       
     
   
   The Hamiltonian is analogous to the Hamiltonian of a mechanical system, which expresses the total energy of a mechanical system. The function y*(τ) is analogous to specifying the lowest energy positions of the components of a mechanical system in a continuous time interval, and the costate function p*(τ) is analogous to specifying the lowest energy momentum of components of a mechanical system in a continuous time interval. Using the above-derived results allows for computation, in a forward direction, along the optimization trajectory, of the partial derivative 
             ∂     y   *         ∂   τ           
The partial derivative
 
             ∂   p       ∂   τ           
can be computed backward from τ=T from the forward-computed y*(τ). Again, with analogy to a classical mechanical problem, the momentum of the components of a mechanical system are essentially determined from the change in position of the components, over time. Thus, momentum, or costate, needs to be computed backwards from forward computed positions.
 
   Minimization of the Hamiltonian is itself an optimization problem, and produces the optimal control policy when the Hamiltonian is minimized. This minimization problem can also be solved using Newton&#39;s approximation method, as follows: 
               ∂   ℋ       ∂   u       =       0   ⁢           ⁢   and   ⁢           ⁢       ∂   ℋ       ∂   r         =       0   ⁢     
     [           u   ⁡     (     τ     k   +   1       )                 r   ⁡     (     τ     k   +   1       )             ]     =       [           u   ⁡     (     τ   k     )                 r   ⁡     (     τ   k     )             ]     -           [               ∂   2     ⁢   ℋ       ∂     u   2                   ∂   2     ⁢   ℋ         ∂   u     ⁢     ∂   r                       ∂   2     ⁢   ℋ         ∂   r     ⁢     ∂   u                   ∂   2     ⁢   ℋ       ∂     r   2               ]       -   1       ⁡     [             ∂   ℋ       ∂   u                   ∂   ℋ       ∂   r             ]       ⁢     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )       ,     r   ⁡     (     τ   k     )       ,     p   ⁡     (     τ   k     )         )                       where                 Φ   κ     ⁡     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )       ,     r   ⁡     (     τ   k     )       ,     p   ⁡     (     τ   k     )         )       =           [               ∂   2     ⁢           ⁢   ℋ       ∂     u   2                   ∂   2     ⁢           ⁢   ℋ         ∂   u     ⁢           ⁢     ∂   r                       ∂   2     ⁢           ⁢   ℋ         ∂   r     ⁢           ⁢     ∂   u                   ∂   2     ⁢           ⁢   ℋ       ∂     r   2               ]       -   1       ⁡     [             ∂   ℋ       ∂   u                   ∂   ℋ       ∂   r             ]       ⁢     (       y   ⁡     (     τ   k     )       ,           ⁢     u   ⁡     (     τ   k     )       ,           ⁢     r   ⁡     (     τ   k     )       ,           ⁢     p   ⁡     (     τ   k     )         )             
Thus, formulation of the minimum time problem and optimal control problem provides a method for computing state vector y and control vectors u and r in a forward direction, and computing the costate vector p in a backward direction, over a range of the continuous iteration variable τ between 0 and a horizon T. This then leads to an iterative optimization method by which the objective function F is forced to a local minimum via an optimal state-vector-trajectory y*(τ) under an optimal control policy u*(τ) and r*(τ).
 
     FIG. 20  graphically illustrates the iterative optimization method based on the above-derived state-vector-change optimization trajectory y* and optimal control policies u* and r*. The iterative optimization method can be thought of as comprising two nested loops, an outer loop controlled by the iteration variable K and the inner loop controlled by the iteration variable k. The outer nested loop represents integration of the partial differentials of the state-vector-trajectory and control policy over a continuous iteration range of 0 to T in order to find a state vector y that minimizes the objective function F. Within each such outer-loop iteration, multiple forward computations of the control policy vectors u and r and the state vector y over small intervals of the continuous iteration variable τ are carried out in steps τ k , where k goes from 0 to a final k value, or, in other words, τ 0 , τ 1 , τ 2 , . . . , τ final , and then, at τ final ≧T, the results of the computations of y(τ k ) and u(τ k ) are used to back calculate the costate vectors p(τ k ) from τ=T to τ=0. During each inner loop, the computation of the y(τ k ) values are based on prior computation of u(τ k ) and r(τ k ), and the computation of both the y(T k ) and u(τ k ) are carried out using costate values p(τ k ) computed at the end of a previous iteration of the outer loop. In the first iteration, the control costate values are approximated as a constant value. 
   In  FIG. 20 , the outer loop is represented by the outer circle  2002 , and the inner-loop computations are represented by step-like radial and circumferential segments, such as step  2004 . The iterative optimization process begins with initial values for the state vector, control vectors, and costate vector of y 0 , u 0 , r 0 , and p 0 , respectively  2006 . The outer-loop iteration variable K is set to 0 ( 2008 ). Then, a number of step-like, inner-loop iterations are carried out, beginning with iteration k=0, each successive iteration of the inner loop incrementing the inner-loop variable k by one. In the first, inner-loop iteration, control vectors U(τ k+1 ) and r(τ k+1)  are computed  2004 . Then, the value of the state vector at τ=τ k+1 , y (τk+1 ), is computed by a number of successive computations  2010 – 2016 . At the end of the first iteration of the inner loop, a determination is made as to whether a local minimum has been reached. If so, then the iterative optimization method terminates, indicated in  FIG. 20  by the dashed line  2018 . Otherwise, a next inner-loop iteration is carried out  2020 – 2021 . Inner-loop iterations continue until τ k ≧T  2022 , at which point the costate vector values p(τ) are computed in a number of steps  2024 – 2030  from τ=T to τ=0, and the outer-loop iteration control variable K is incremented to conclude a full outer-loop iteration and return to the beginning point  2032  for a next outer-loop iteration. 
   Considering  FIGS. 19 and 20  together provides a graphical, conceptual view of the iterative optimization process. Each inner-loop iteration of the optimization process can be considered to correspond to a segment in the optimization trajectory shown in  FIG. 19 , such as segment  1906 . One full outer-loop iteration, from τ=0 to τ=T can be considered to produce a large segment of the optimization trajectory, such as the optimization trajectory segment beginning at initial point  1904  and extending down to point  1908 . Multiple outer-loop iterations therefore compute a number of large optimization trajectory segments that eventually drive the state vector y to a value at which the objective function F is at least a local minimum.  FIG. 21  provides additional, graphical illustration of the overall iterative optimization method. As shown in  FIG. 21 , a first outer-loop iteration produces the optimization trajectory segment from an initial point  2102  to point  2104 . A second outer-loop iteration computes the optimization trajectory segment  2106  from point  2104  to point  2108 . Third outer-loop iteration produces the optimization trajectory segment  2110  from point  2108  to local optimum  2112 . Note that the value of the horizon, T, may change with each outer-loop iteration. Thus, in  FIG. 21 , the endpoint of each optimization trajectory segment is marked T K , with K=1, 2, and 3, indicating that the endpoint of each outer-loop iteration serves as the initial point for each subsequent outer-loop iteration. It should also be noted that the magnitude of Δτ for each inner-loop step may also vary. 
     FIG. 22  is a control-flow diagram that describes, in general, the iterative optimization step ( 203  in  FIG. 2 ). First, initial values for y, u, r, and p are chosen, the outer-loop iteration variable K is initially set to 0, and an initial horizon T 0  is computed in step  2202 . Next, in an inner loop comprising steps  2204 – 2208 , the control vectors u and r and state vector y are propagated forward from τ=0 to τ=T K  in τ increments of τ k . First, in step  2204 , a next value for the control vectors u and r is computed. Then, in step  2205 , the state vector y is propagated from τ=τ k  to τ=τ k +Δτ k  in step  2205 . In step  2206 , a determination is made as to whether the optimization has converged. If so, then the iterative optimization has completed, and returns. Otherwise, a determination is made as to whether τ k +Δτ remains less than T K , in step  2207 . If not, then τ k+1  is assigned to τ k +Δτ, and the inner-loop control variable k is incremented, in step  2208 . However, if τ k  is greater than or equal to T K , then the costate values are computed backwards from T K  to 0, in step  2209 , and the iteration variables and a database containing the computed values of the state vector, control vectors are updated in step  2210 . 
   In convex problems, in which there is a single critical point, the iterative optimization method described above with reference to  FIGS. 19–21  normally can converge in a single thread within a reasonable amount of time even for extremely large hyper dimensional problem domains with thousands of decision variables. However, optimization problems in applications are usually not convex, and the problem domains may feature many critical points.  FIG. 23  illustrates one such problem domain. In  FIG. 23 , the surface boundary for the problem domain  2303  is bowl-shaped with a number of prominent central protrusions, such as protrusion  2304 . These protrusions may be thought of as a small mountain range with three peaks ( 2304 ,  2306  and  2308 ) and a deep local valley, or bowl  2310 . The global optimization amount for the problem domain lie deeply below the mountain range in the depths of the outer bowl. 
     FIG. 24  shows the same problem domain as shown in  FIG. 23 , or critical points explicitly indicated. The peaks of the mountain range  2304 ,  2306  and  2308  represent local maxima within the problem domain, and the low point in the small valley  2310  represents a local minimum. Saddle points within the mountain range, such as saddle points  2312  and  2314 , also represent critical points.  FIG. 25  shows a magnified view of a critical point, such as critical point  2304  in  FIG. 24 . As shown in  FIG. 25 , the partial differentials calculated for a number of points  2502 – 2504  near the critical point in the x and y directions calculated for a number of points  2502 – 2504  near the crucial point  2506  are all small in absolute value. Therefore, the gradient calculated for the points  2202 – 2304  near the critical point  2506  are also close to 0. To recapitulate, the derivative of the objective function with respect to τ, 
                 (       ∂     F   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )       ,     r   ⁡     (   τ   )         )           ∂   y       )     T     =   G     ,         
gives a means for determining the optimization trajectory near a critical point. Therefore, when the gradient is 0 or near 0, as it is for points near a critical point, then the rate of the state-vector trajectory with respect to τ is also 0. Moreover, G is sensitive near the critical point, since small local variations can cause relatively large magnitude changes in G and can even switch the sign of G. However, when G is small, then Newton&#39;s method, and other similar methods used to compute state-vector trajectories, are generally not useful for finding an optimization trajectory. In other words, the gradient field generally provides the direction for near-optimal or optimal trajectories, and therefore, while the gradient field is 0, all possible trajectories appear to be equally good. It is a situation analogous to attempting to use a compass at the magnetic north or south poles. As will be discussed later, a different method for computing the optimization trajectory for the state vector needs to be employed near critical points.
 
   There is another consideration with regard to critical points. This consideration involves search-like strategy that needs to be employed near critical points by a global optimization method in order to attempt to identify and evaluate a number of potential minima within a problem domain.  FIG. 26  illustrates a first-search technique. In  FIG. 26 , an iteration of the outer-loop of the optimization method to horizon T n  has carried the state vector to point  2602 . This point is near critical point  2306 , the peak of the mountain range. The fact that the gradient is nearly 0 at this point informs the optimization method that it is arrived at the vicinity of a critical point. However, its next step depends on the type of critical point at which it has arrived. In the case of a local maximum, such as point  2306 , the optimization method needs to employ some method for perturbing the optimization trajectory away from the local maximum  2306  to continue a descent  2308  towards a local minimum. However, in the neighborhood of a local minimum, such as local  2310 , the optimization method needs to employ a technique to ascend out of the local minimum  2312  to another point  2314  from which the optimization method can continue to descend  2316  towards a local minimum. 
     FIG. 27  illustrates a second strategy for searching a problem domain. In  FIG. 27 , the search method has again arrived at a point  2602  near a local maximum  2306 . However, upon detecting the local maximum, the optimization method spawns a number of additional, relatively independent optimization threads which continue descending in different directions  2702 – 2705 . One of the relatively independent threads  2704  reaches a point  2708  near a saddle point  2314 . At this point, this thread,  2704  forks, or spawns, an additional thread, so that each of the resulting relatively independent threads descends from the saddle point in different directions  2710  and  2712 . Another of the originally spawned, relatively independent threads  2702  reaches a local minimum  2310 . At this point, rather than spawning new optimization threads, the relatively independent thread  2702  terminates, reporting the local minimum back to its parent thread  2701 . Parent threads can select the lowest local minimum identified by their children and return that local minimum to their parents, resulting in a recursive search of the problem domain to identify a near optimal or optimal solution to the optimization problem. In alternative embodiments, parent threads can use information about reported local minima to truncate the problem domain, in order to modify the problem domain based on the local minima, and communicate the truncation of the problem domain to all other relatively independent threads. 
   By either the single-threaded technique shown in  FIG. 26 , or the multi-threaded technique shown in  FIG. 27 , the optimization method needs to have a means of determining the type of critical point that is arrived at when the gradient falls to a value near 0. The critical-point beacon for illuminating types of critical points is the Hessian matrix, described above. When the gradient is near 0, then when all elements of the Hessian matrix are negative, a local maximum has been reached. When the gradient is 0 and all elements of the Hessian matrix are positive, then a local minimum has been reached. Otherwise, when the gradient is 0, and the elements of the Hessian matrix have mixed signs, then an inflection point or saddle point has been reached. The critical-point beacon properties of the non-singular Hessian matrix are formerly described as follows: 
   for a point y, if G| y ≈0 
   
       
       
         
           if all eigenvalues of H are negative, then y is near a local or global maximum 
           if all eigenvalues of H are positive, then y is near a local or global minimum otherwise y is near a saddle point 
         
       
     
  
   A single optimization thread, as discussed above with reference to  FIGS. 26 and 27 , comprises a collection of hierarchically organized modules that may execute concurrently in separate computational environments in order to distribute the large computational tasks involved in propagating the state vector according to the above-described optimization method.  FIG. 28  illustrates the modules and modular organization of an optimization thread. The highest-level module is referred to, in this document, as a “super”  2802 . The super includes a relatively large database  2804  of intermediate results. The super is responsible for computing each successive horizon T K , for choosing an initial point, in the case of the root thread in a multi-threaded optimization method, as described above with reference to  FIG. 27 , and for launching and coordinating the activities of subordinate modules. The first-level subordinate modules include a u_and_r module  2806 , a y module  2808 , and a p module  2810 . The u_and_r module  2806  is responsible for the forward computation of the control vectors u_and_r, as shown in step  2204  in  FIG. 22 . 
   The y module  2808  is responsible for propagation of the state vector through successive intervals Δτ, shown as step  2205  in  FIG. 22 . The p module  2810  is responsible for computing costate vector values backwards from T k  to 0 at the end of each outer-loop iteration of the iterative optimization method, shown as step  2209  in  FIG. 22 . 
   The u_and_r module  2806  distributes the control vector forward calculation amongst m agents  2812 – 2815 . They module  2808  and the p module  2810  each distributes its respective task amongst n agents  2816 – 2818  and  2820 – 2823 , respectively. Distribution of the y propagation by y module  2808  over the n y-module agents  2816 – 2818  is an important approach to decomposing the overall optimization computations and providing opportunities for parallel computing. The state-vector propagation computations undertaken by the y module  2808  and its respective agents  2816 - 2818  are the central, and most computationally expensive, tasks involved in computing the state-vector-change optimization trajectory. In an alternative embodiment of the present invention, the p module uses an efficient, interpolation method for estimating costate vector values, and thus does not need to distribute the problem among n p-module agents. Thus, the p-module computes costate without distribution. In other alternative embodiments, the p-module may distribute costate interpolation over a small number of p-module agents. 
   The database  2804  maintained by the super  2802  is large, and may be contained not only in memory, but in mass storage devises directly or remotely coupled to one or more computer systems over which the modules of a thread are distributed.  FIG. 29  illustrates, in part, the contents of the database maintained by the super in various embodiments of the present invention. In general, the super must maintain the values of various quantities calculated by the modules at each subinterval rk of the horizon from 0 to T. Thus, in  FIG. 29 , each sub-column of the first column of data structures, such as sub-column  2902  of the vector G data structure  2904 , represents the value of the corresponding computed entity at a particular value of τ. In  FIG. 29 , sub-column  2902  represents the value of the vector G at τ=0. The super stores values of the state vector  2906 , the control vectors u and r  2908 – 2909 , the vector w  2910 , the scaler ρ  2911 , a Morse-flags vector  2912 , the costate vector p  2914 , the scaler values of the objective function F  2915 , the vector G  2904 , the value of the Hamiltonian  2916 , the current T and iteration control variables and other such local control variables  2918 , and the inverses of the Hessian matrix  2920 – 2927 . 
     FIGS. 30–34  illustrate, using the illustration conventions of  FIG. 28 , the flow of information among the hierarchically organized modules of an optimization thread during one outer-loop iteration of the optimization method. In the case of a root thread, as shown in  FIG. 30 , an initial point for the optimization trajectory is determined  3002 . In addition, initial values for the state vector, costate vector, control vector, and control variables, among other data structures, are determined for initializing the various data structures and elements  3004 . Non-root threads inherit a starting point and current values for the various computational entities, such as the state and control vectors, from parent threads. As shown in  FIG. 31 , in the first step of the inner loop ( 2204  in  FIG. 22 ), the u_and_r module  2806  receives the values of the state vector, costate vector, vector w, and the inverse of the Hessian from the super  3102  and distributes these values and/or portions of these values to u_and_r-module agents  3104 – 3107 , receives intermediate computed results back from the u_and_r-module agents  3108 ,  3109 ,  3110 , and  3111 , computes a new control vector value, and returns the new control vector value  3112  back to the super.  FIG. 32  illustrates data flow during the next step of the inner loop ( 2205  in  FIG. 22 ). In this step, they module  2808  receives the newly computed control vectors u and r from the super  3202 , distributes these values, along with additional values obtained from the super and/or computed by the y module, including the Hessian, vector other such values, to the y agents  3204 – 3206 , receives intermediate computed values back from the y agents  3208   3209   3210 , and computes and returns to the super computed values for the state vector, w vector, differential of the w vector with respect to the state vector, and the inverse of the Hessian, among others  3210 . The y module carries out a number of individual propagation steps for each inner-loop iteration, as discussed above.  FIG. 33  indicates that the data-exchange steps illustrated in  FIGS. 31 and 32  are successively repeated in a number of inner-loop iterations.  FIG. 34  illustrates the data exchange involved in costate vector calculations undertaken by the p module  2810  in step  2209  of  FIG. 22 . The p module receives computed values of the state vector, control vectors, and other such values from the super  3402 , distributes these values, or portions of the values, among the p agents  3404 - 3407 , receives intermediate results from the p agents  3408 – 3411 , and returns computed values for the costate vector over a range of τ values from 0 to T K    3412  to the super. As discussed above, in currently preferred embodiments, the p module interpolates, rather than integrates, costate vectors, and thus does not need to distribute the problem in the manner illustrated in  FIG. 34 . 
   In the following discussions, for notational convenience, the u and r vectors are considered together as the u vector. In other words, in the subsequent discussion, the binding variable coefficients and barrier-function coefficients together compose a single control vector u. Recall that the iterative optimization strategy employed by embodiments of the present invention involves integrating 
                 ⅆ   y       ⅆ   τ       ⁢     (   τ   )       ,         
expressed as follows:
 
                 ⅆ   y       ⅆ   τ       ⁢     (   τ   )       =     -       w   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )         ρ   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )                     where               w   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )       =         (     H   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )       )       -   1       ⁢     G   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )                       ρ   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )       =         (     G   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )       )     T     ⁢     w   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )               
In one embodiment, the integration of
 
               ⅆ   y       ⅆ   τ       ⁢     (   τ   )           
is carried out by a forward Euler integration method with variable step size Δτ k :
 
               y   i     ⁡     (       τ   k     +     Δτ   k       )       =         y   i     ⁡     (     τ   k     )       +     Δ   ⁢           ⁢     τ   k     ⁢       ⅆ     y   i         ⅆ   τ       ⁢     (     τ   k     )               
Therefore, as is apparent in the above expressions, they module, for each forward integration step, needs to invert the Hessian matrix in order to compute the w vector as the product of the inverse of the Hessian and the vector G  FIG. 35  graphically illustrates this computation. Inversion of the Hessian matrix is by far the most computationally intensive task involved in optimization. Therefore, the y module distributes inversion of the Hessian matrix among the y-module agents, relying on several assumptions for this distribution.
 
   Recall that: 
             H   ⁡     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )       ,     w   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )         )       =         [       ∂   G       ∂   y       ]     T     ⁢     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )         )               FIG. 36  illustrates multiplication of the Hessian matrix by the w vector to produce the g vector in a distributed fashion. As shown in  FIG. 36 , the G and w vectors, and the Hessian matrix, can be vertically partitioned into partitions i=0 to 5. The first partition of the G vector  3602  is referred to as G 0 , the first vertical partition of the Hessian matrix  3604  is referred to as H 0 , and the first partition of the w vector  3606  is referred to as w 0 . The Hessian matrix can additionally be partitioned horizontally into horizontal partitions j=0 to 5. In  FIG. 36 , the first horizontal partition of the Hessian matrix  3608  appears as a first vertical column within the Hessian matrix.
 
   Both vertical and horizontal partitioning of the Hessian matrix results in a block-like partitioning of the Hessian matrix into submatrices. For example,  FIG. 36  shows the Hessian matrix partitioned into thirty-six submatrices, beginning with the left-most and upper-most submatrix H 0,0    3610  and ending with the right-most and lower-most submatrix H 5,5    3612 . Note that, in the case of a diagonally dominant Hessian matrix, in which the diagonal block submatrices H 0,0 , H 1,1 , H 2,2 , H 3,3 , H 4,4 , and H 5,5  have elements with magnitude much greater than the elements of other non-diagonal submatrices, then matrix multiplication of the Hessian matrix by the w vector can be approximated by matrix multiplication of w-vector partitions by the block-diagonal submatrix within the corresponding vertical partition of the Hessian matrix. For example, in a full matrix multiplication, the G-vector partition G 0  equals the sum of the matrix-multiplication products of the partitions of the w vector with corresponding blocks of the first vertical partition of the Hessian matrix. Formally, 
   
     
       
         
           
             
               
                 
                   G 
                   i 
                 
                 = 
                 
                   
                     H 
                     i 
                   
                   ⁢ 
                   w 
                 
               
             
           
           
             
               
                 = 
                 
                   
                     ∑ 
                     
                       j 
                       = 
                       0 
                     
                     
                       J 
                       - 
                       1 
                     
                   
                   ⁢ 
                   
                     
                       H 
                       ij 
                     
                     ⁢ 
                     
                       w 
                       j 
                     
                   
                 
               
             
           
           
             
               
                 
                   H 
                   i 
                 
                 = 
                 
                   [ 
                   
                     
                       
                         
                           H 
                           
                             i 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             0 
                           
                         
                       
                       
                         ⋯ 
                       
                       
                         
                           H 
                           
                             i 
                             ⁡ 
                             
                               ( 
                               
                                 J 
                                 - 
                                 1 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                   ] 
                 
               
             
           
         
       
     
   
   When the Hessian matrix is diagonally dominant, then the off-diagonal submatrices have relatively less magnitude and are of less importance in computation of a w-vector partition. This situation occurs when the Hessian matrix is partitioned in a way that removes, as much as possible, dependence between each partition and all other partitions. Therefore, a particular y-module agent can approximately calculate the current value of a w-vector partition  3810  by inverting a corresponding Hessian block-diagonal submatrix  3808  and multiplying the large, bracketed term  3812  in  FIG. 38  by the inverse of the block-diagonal Hessian submatrix  3808 , using previously calculated values for the remaining w-vector partitions  3814 – 3817 . 
     FIG. 39  illustrates computation of a w-vector partition w 1  by the second y-module agent. The second y-module agent  3902  obtains the current Hessian matrix block-diagonal submatrix H 1,1  and current G-vector partition G 1  and most recently calculated values of the w-vector partitions w 0 , w 2 , w 3 , w 4 , and w 5  from the y-module  2808 , computes the inverse of the current Hessian block-diagonal submatrix H 1,1  and then, using the computed inverse of the Hessian matrix, computes the current approximate value of the w-vector partition w 1 . The y-module agent  3902  can then, in a next computational step, fumish the approximate current value of the w-vector partition to the y-module  2808  for distribution to the other y-module agents, and to receive their most recent computed values for the other w-vector partitions w 0 , w 2 , w 3 , w 4 , and w 5  in order to compute a next w-vector partition value w 1  at the next τ increment. Using this method, the y-module agents can synchronously compute the approximate value of w-vector partitions at each step. 
   Next, a control-flow diagram illustration of the super module is provided.  FIG. 40  is a control-flow diagram illustrating a highest conceptual level of the super module. In step  4002 , the super module determines whether or not it is a replicate, or child thread, spawned or forked from a parent thread. If not, and the current thread is therefore the root thread, then in step  4004 , the super module determines an initial point for the optimization and, in step  4006 , initializes various computational entities, including the state vector, costate vector, and other sets of values, and also initializes a number of variables, including the Boolean variables “critical” and “precritical” to the Boolean value FALSE. Following initialization by a root thread, or as a result of flow of control from step  4002  for a replicate thread, the super module next chooses, in step  4008 , a horizon T. Then, the super module invokes the function “integrate_u_and_y_to_T,” in step  4010 , which carries out successive inner-loop interactions of the optimization method, described above. Next, in step  4012 , the super module invokes the routine “evaluate” in order to make decisions based on the inner-loop iterations carried out by the routine “integrate_u_and_j_to_T” in step  4010 . The routine “evaluate” returns a Boolean value “continue” which the super module evaluates, in step  4014 , in order to determine whether or not to carry out yet another iteration of the outer loop of the optimization method, described above. If so, then the super module calls the routine “integrate p”  4016  to back calculate costate-vector values and then returns to step  4008  to initiate another outer-loop iteration. If the super module decides not to continue outer-loop durations, then the super module invokes the routing “report”, in step  4018  and then returns. 
   Control-flow diagrams illustrating steps  4004 ,  4008 ,  4010 ,  4012 , and  4018  are next provided. The routine “integrate p” called in step  4016 , and the routines “integrate u” and “integrate y,” called from the routine “integrate_u_and_to_T” are subsequently described mathematically, rather than using control-flow diagrams, in the interest of clarity and brevity. 
     FIG. 41  is a control-flow diagram for the routine “initial point” called in step  4004  of  FIG. 40 . In order to determine an initial value for the state vector, y 0 , the initial point routine, in step  4102 , guesses at a value using any supplied hints or heuristics from users or from the mathematical modeling or preprocessing portions of the optimizer. Then, in the for-loop comprising steps  4104 – 4107 , the initial point routine evaluates each constraint g i  in order to determine whether or not the constraint, evaluated at the initially guessed state-vector point y 0 , is greater than 0. If so, then the constraint is placed into partition I in step  4107 . If not, then the constraint is placed in partition J in step  4106 . Then, in step  4110 , the initial point routine constructs an objective function {tilde over (F)} as follows: 
             F   ~     =       -       ∑     i   ∈   I       ⁢       u   i     ⁢     log   ⁡     (       g   i     ⁡     (   y   )       )             -       ∑     i   ∈   J       ⁢       u   j     ⁢     log   ⁡     (         g   j     ⁡     (   y   )       +     T   j       )           -       ∑     l   ∈   J       ⁢       u   l     ⁢     log   ⁡     (     T   l     )                   
which is then minimized in step  4112  by a recursive call to the optimizer ( FIG. 30 ) to produce a modified initial state-vector value. The objective function {tilde over (F)} can be thought of as initially expanding the problem-domain boundary of the objective function F in order to envelope the initial value for the state vector, and the optimization of {tilde over (F)} then contracts the problem-domain boundary, forcing the initial point into the interior of the problem domain of the objective function F. In step  4114 , the initial point routine determines whether the modified state-vector value is feasible. A feasible state-vector value is an interior point of the problem domain with a reasonably high probability of being close enough to a local minimum for the optimization method to find a local minimum. If the modified state-vector value is feasible, then it is returned as the initial point for the optimization. Otherwise, control flows back to step  4102  for another pass at determining an initial state-vector value. In subsequent iterations, information gleaned from previous iterations may be employed in step  4102  in order to better guess an initial value for the state vector.
 
     FIG. 42  is a control-flow diagram for the routine “choose_next_T.” In step  4201 , the routine “choose_next_T” assigns local variable “T” to a minimum horizon constant. Then, in the for-loop, comprising steps  4203 – 4207 , the routine “choose_next_T” calculates, for each y-module agent, a T value obtained by dividing a constant parameter by the minimum absolute eigenvalue of the block-diagonal Hessian submatrix associated with the y-module agent. If the computed T value for the y-module agent is greater than the contents of the variable “T,” as detected in step  4205 , then the variable “T” is assigned to the T value calculated for they-module agent. Next, in step  4208 , the routine “choose_next_T” initializes local variables “n” to a maximum number of increments, initializes local variable “a” to an increment coefficient, and initializes a loop-control variable “j” to  0 . Then, in steps  4210 – 4212 , the routine “choose_next_T” increases the value of the variable “T” in step  4212  until the derivative of the costate vector with respect to T evaluated for the current contents of the variable “T” approaches 0, as detected in step  4210 . The loop may also terminate when the number of iterations is equal to the contents of the local variable “n.” When 
                 ⅆ   p       ⅆ   T       =   0     ,         
then computation of the costate vectors over the continuous iteration range 0 to T can be more easily calculated. In alternative embodiments of the present invention, the costate vectors can be computed using a polynomial approximation, rather than employing complex integration.
 
     FIG. 43  is a control-flow diagram of the routine “evaluate,” called in step  4012  of  FIG. 40 . In step  4302 , the routine of “evaluate” sets the Boolean flags “at_minimum” and “continue” to the Boolean value FALSE. Next, in step  4304 , the routine evaluate determines whether the computational resources allocated for the optimization have been exhausted. This may include continuous computation for longer than a maximum specified time, the usage of processing resources, memory, or a combination of computational resources greater than a specified maximum amount of resources to be used, the occurrence of error conditions, or other types of predetermined termination criteria. If the computational resources are exhausted, then the routine evaluate returns. Otherwise, the routine evaluate, in step  4306 , determines whether the optimization thread is currently critical. An optimization thread is critical if all y-module agents are currently in Morse mode, and have been so for some predetermined number of iterations. Morse mode is described in detail, below. If the current thread is not critical, then the routine “evaluate” sets the Boolean flag “continue” to TRUE, in step  4308 , and then returns. Otherwise, in step  4310 , the routine “evaluate” determines whether the thread, currently determined to be critical, has reached a local minimum. If so, then the Boolean flag “at-minimum” is set to TRUE, in step  4312 , and the routine “evaluate” returns. If the optimization thread is not at a local minimum, but is critical, then the routine “evaluate,” in step  4314 , determines a number of replicant threads to spawn, and also initial directions or starting points for the replicant threads. Spawning of replicant threads is illustrated in  FIG. 27 . In the for-loop of steps  4316 – 4319 , the determined number of replicant threads are spawned, each replicant thread linked has a child to the current optimization thread. Then, in step  4320 , the routine “evaluate” may attempt to perturb the current optimization trajectory in order to move the optimization trajectory away from the critical point and continue descending towards a local minimum. In step  4322 , the routine “evaluate” sets the Boolean flag “continue” to the Boolean value TRUE that then returns. 
     FIG. 44  is a control-flow diagram for the routine “report.” This routine is invoked in step  4018  of  FIG. 40 . In step  4402 , the routine “report” determines whether the current optimization thread has a parent. If not, and the current optimization thread is therefore the root optimization thread, then the routine “report” determines, in step  4404 , whether there are any children threads still running that are linked to the current, root optimization thread. If so, then the current, root optimization thread waits, in step  4406 , for the children optimization threads to finish execution. Next, in step  4408 , the routine “report” sets local variable “min” to the current value of the objective function computed with the current state-vector value, sets the local variable “location” to the current state-vector value, and sets the local variable “at min” to the Boolean value flag “at_minimum” returned by the routine “evaluate” for the current optimization thread. In step  4410 , the routine “report” determines whether any children threads have terminated and reported back to the current optimization thread. If so, then in the for-loop of steps  4412 - 4415 , the current optimization thread evaluates the report generated by each terminated child in order to determine the lowest value for the objective function reported by any child and achieved by the currently executing optimization thread. Next, in step  4416 , the routine “report” determines whether there are any children still running. If so, then in thefor-loop of steps  4418 – 4420 , the routine “report” links each child to the parent of the currently executing thread. Note that the test in step  4416  fails in the case that the currently executing optimization thread is the root thread. If the currently executing optimization thread is the root thread, as determined in step  4422 , then the routine “report” returns the minimum objective function value, state vector, and indication of whether the optimization is currently at a minimum in step  4424 . Otherwise, the routine “report” reports these same values back to its parent thread, in step  4426 . 
     FIG. 45  is a control-flow diagram for the routine “integrate_u_and_to_T,” called in step  4010  of  FIG. 40 . In the for-loop of steps  4502 – 4507 , the routine “integrate_u_and_y_to_T” successively calls the routines “integrate u” and &#39;steps  4503  and  4504 , in order to complete one inner-loop execution for the optimization method. The successive iterations continue until either the optimization thread goes critical, as detected in step  4505 , until the value τ k  is greater than or equal to T, or until some other stopping criteria are detected, such as exhaustion of computational resources or other such stopping criteria in step  4507 . 
     FIG. 46  is a control-flow diagram for the routine “integrate y.” This routine describes the operation of the y module ( 2808  in  FIG. 28 ). In step  4602 , the routine “integrate y” obtains symbolic representations of the Hessian matrix, the vector G, and other data needed by the routine “integrate y” from the super module. Next, in step  4604 , the routine “integrate y” computes numerical values for the Hessian matrix and vector G. In step  4606 , the routine “integrate y” distributes partitions of the computed Hessian and G vector to each y-module agent. The y-module agents compute their associated w-vector partition values and distributed row values, which are collected by the routine “integrate y” in step  4608 . In step  4610 , the routine “integrate y” determines whether any agents are currently in Morse mode. 
   The term “Morse mode” refers to an agent that currently needs to compute its w-vector partition based on the Morse quadratic approximation of the original objective function associated with the y agent method since the agent is currently near a critical point with respect to the agent&#39;s partition of the vector G. If any agents are in Morse mode, then, in step  4612 , the routine “integrate y” directs the Morse-mode agents i to recompute p i  and w i  and collects the recomputed p i  and w i  from each Morse-mode agent i. Next, in step  4614 , the routine “integrate y” determines whether all agents are currently in Morse mode. If so, then, in step  4616 , the routine “integrate y” determines whether the local variable “precritical” currently has the Boolean value TRUE. If so, then the agents were previously in Morse mode in the previously executed iteration, and so the local variable “critical_count” is incremented, in step  4618 . In step  4619 , the routine “integrate y” compares the contents of the local variable “critical_count” to a threshold constant and, when the contents of the local variable “critical_count” equals or exceeds the critical constant, the routine “integrate y” concludes that the current optimization thread has remained critical for a sufficient number of iterations to declare the optimization thread to be critical by setting the Boolean flag “critical” to Boolean value TRUE in step  4620 . However, if the contents of the variable “critical_count” do not justify declaring the current optimization thread to be critical, then control flows to step  4624 . If, back in step  4616 , the routine “integrate y” determines that the variable “precritical” did not have the Boolean value TRUE, then, in step  4622 , the variable “precritical” is set to true and the variable “critical_count” is initialized to 0 before control flows to step  4624 . In step  4624 , the routine “integrate y” computes global ρ and w-vector values, and it computes from them the value of 
   
     
       
         
           
             
               ⅆ 
               y 
             
             
               ⅆ 
               τ 
             
           
           . 
         
       
     
   
   Next, in step  4625 , the routine “integrate y” determines a value Δτ for the current Euler forward integration step. If the computation of Δτ succeeds, as determined in step  4626 , then the routine “integrate y” carries out the forward integration step, in step  4628 , updates the binding variable coefficients, if necessary, in step  4630 , and returns. If the computation of Δτ fails, as detected in step  4626 , then the routine “integrate y” returns with an error. 
   The forward Euler integration technique has been discusses above, but is presented formally below. The global value for ρ is assembled from the partitioned ρ values computed by y-module agents in a current inner-loop step as follows: 
             ρ   ⁡     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )         )       =       ∑     i   =   0       n   -   1       ⁢       ρ   i     ⁡     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )         )               
Next, the values for the partitions of the w vector are collected from the y-module agents and assembled to produce a w vector. Then, a continuous iteration increment Δτ is computed, by a method to be described below. Finally, the state vector is propagated from a current continuous iteration value τ k  to a next continuous iteration value τ k+Δτ  as follows:
   y   i (τ k   +Δτ   k )= y   i (τ k )+Δ k   dy   i   l /d τ(τ k ) 
Finally, the continuous iteration value τ k  is updated as follows:
 τ k =τ k +Δτ k   
Note that ay-module agent is considered to be critical with respect to the y-module&#39;s vector-G partition when the absolute value of the y-module agent&#39;s local ρ, ρ i , falls below some constant parameter. The computation of a w-vector partition w i  and the y-module-agent-local ρ, ρ i , is described below.
 
     FIG. 47  is a control-flow diagram for the routine “compute Δτ,” called in step  4625  of  FIG. 46 . In step  4702 , a guess, or estimate, Δτ′, is made for Δτ. This guess or estimate may be obtained from a constant parameter, or may be, in more elaborate embodiments, derived by empirical methods or heuristics from prior computed values of Δτ and from recent history of state-vector values, control-vector values, and other such information. In the described embodiment, for forward Euler integration steps, when the control vector has an approximate constant value, then F decreases by one at each step: 
                 ⅆ   F       ⅆ   τ       ⁢     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )         )       =     -   1.           
When some of the y-module agents are in Morse mode,
 
             ⅆ   F       ⅆ   τ           
may depart from the value “−1.” However, if the value of
 
             ⅆ   F       ⅆ   τ           
ever rises above 0, as detected in step  4704  of the routine “compute Δτ,” then an error condition obtains since
 
               ⅆ   y       ⅆ   τ       ⁢     (     τ   k     )           
is not a descent direction and therefore the forward Euler integration method cannot produce a feasible trajectory y*(τ). Otherwise, the estimate of Δτ, Δτ′, is inserted into the forward Euler integration step equation
 
             y   ⁡     (       τ   k     +     Δ   ⁢           ⁢     τ   ′         )       =       y   ⁡     (     τ   k     )       +     Δ   ⁢           ⁢     τ   ′     ⁢       ⅆ   y       ⅆ   τ       ⁢     (     τ   k     )               
and the state vector is propagated to y(τ k+  Δτ′). This computed next state-vector value is checked, in step  4706 , for feasibility. If the state-vector value is not feasible, then the estimate of Δτ, Δτ′, is decreased, in step  4708 , and a new forward state-vector value is again computed and checked for feasibility in step  4706 . Once a feasible forward state-vector is found with a corresponding Δτ estimate, Δτ′, then, in step  4710 , the Armijo condition
 
               F   ⁡     (       y   ⁡     (       τ   k     +     Δ   ⁢           ⁢     τ   ′         )       ,     u   ⁡     (     τ   k     )         )       -     F   ⁡     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )         )         &lt;     σ   ⁢           ⁢   Δ   ⁢           ⁢     τ   ′     ⁢       G   T     ⁡     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )         )       ⁢       ⅆ   y       ⅆ   τ       ⁢     (     τ   k     )             
is evaluated with respect to the estimate of Δτ, Δτ′. If the Armnijo condition is satisfied, then the estimate of Δτ, Δτ′, is returned as the new computed Δτ. Otherwise, Δτ is approximated. Several different approximation methods are available. In the case that a previous feasible estimate of Δτ, Δτ A , is available, then a quadratic approximation to the objective function F around τ k  can be made, as follows:
   F ( y (τ k +Δτ),  u (τ k ))≈ {tilde over (F)} (Δτ)= aΔτ   2   +bΔτ+c   
Then, an estimate for Δτ is obtained as follows:
 
               F   ~     ⁡     (   0   )       =       F   ⁡     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )         )       =   c                     F   ~     ⁡     (     Δ   ⁢           ⁢     τ   A       )       =       F   ⁡     (       y   ⁡     (       τ   k     +     Δ   ⁢           ⁢     τ   A         )       ,     u   ⁡     (     τ   k     )         )       =         a   ⁡     (     Δ   ⁢           ⁢     τ   A       )       2     +     b   ⁡     (     Δ   ⁢           ⁢     τ   A       )       +   c                         ⅆ     F   ~           ⅆ   Δ     ⁢           ⁢   τ       ⁢     (   0   )       =           (     G   ⁡     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )         )       )     T     ⁢       ⅆ   y       ⅆ   τ       ⁢     (     τ   k     )       =   b                 c   =       F   ~     ⁡     (   0   )                   b   =         (     G   ⁡     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )         )       )     T     ⁢       ⅆ   y       ⅆ   τ       ⁢     (     τ   k     )                   a   =           F   ~     ⁡     (     Δ   ⁢           ⁢     τ   A       )       -       F   ~     ⁡     (   0   )       -     Δ   ⁢           ⁢         τ   A     ⁡     (     G   ⁡     (       y   ⁡     (     τ   k     )       ,     u   ⁡     (     τ   k     )         )       )       T     ⁢       ⅆ   y       ⅆ   τ       ⁢     (     τ   k     )             (     Δ   ⁢           ⁢     τ   A       )     2             
A Δτ along the descent direction that minimizes F can be approximated by finding a Δτ that minimizes {tilde over (F)}, which occurs when:
 
               ⅆ       F   ~     ⁡     (     Δ   ⁢           ⁢   τ     )           ⅆ     (     Δ   ⁢           ⁢   τ     )         =         2   ⁢   a   ⁢           ⁢   Δ   ⁢           ⁢   τ     +   b     =   0           
Therefore:
 
             Δ   ⁢           ⁢     τ   ′       =     -     b     2   ⁢   a               
When two previous feasible estimates for Δτ for the current τ k  are known, then a cubic approximation method may be employed. However Δτ is approximated in step  4712 , the new approximate Δτ is used to recompute a forward value for the state vector, which is then tested for feasibility in step  4706 .
 
   As discussed above, when a y-module agent is near a critical point and is therefore in Morse mode, the standard forward Euler integration procedure discussed above is not useful in finding an optimal optimization trajectory y*. In general in the neighborhood of a critical point, the Hessian block-diagonal submatrix H i,i  corresponding to ay-module agent i may have negative eigenvalues, and therefore the computed direction might not be a descent direction. However, the Hessian submatrix H i,i  can be modified in such a way as to produce a positive definite Hessian submatrix Ĥ i,i  in order to carry out a Morse-mode forward integration of the state vector. 
   The derivation of the Morse-mode approach to forward integration of the state vector is based on the following three assumptions: (1) the objective function F(y(τ), u(τ)) is roughly constant with respect to u(τ); (2) the Hessian matrix H(y(τ), u(τ)) is diagonally dominant, implying that the y-module agents are partitioned for relative independence; and (3) the state vector y i (τ) changes very little while y-module agent i is in Morse mode. Assume that y-module agent i enters Morse mode, or local criticality, at τ k ={circumflex over (τ)} i , as determined by the value of ρ i (y(τ), u(τ)) calculated by the standard method. Then according to the first and third assumptions, above, the objective function F may be approximated about the domain point (y 0 (τ k ), . . . , y 0 ({circumflex over (τ)} i ), . . . y n−1 (τ k ), u({circumflex over (τ)} i )) by a Taylor approximation as follows: 
             F   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )       =       F   ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       +         (       G   i     ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       )     T     ⁢     (         y   i     ⁡     (   τ   )       -       y   i     ⁡     (       τ   ^     i     )         )       +       ∑       l   =   0       l   ≠   i         n   -   1       ⁢         (       G   l     ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       )     T     ⁢     (         y   l     ⁡     (   τ   )       -       y   l     ⁡     (       τ   ^     i     )         )         +       1   2     ⁢       (         y   i     ⁡     (   τ   )       -       y   i     ⁡     (       τ   ^     i     )         )     T     ⁢       H   ii     ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       ⁢     (         y   i     ⁡     (   τ   )       -       y   i     ⁡     (       τ   ^     i     )         )       +         (         y   i     ⁡     (   τ   )       -       y   i     ⁡     (       τ   ^     i     )         )     T     ⁢       ∑       l   =   0       l   ≠   i         n   -   1       ⁢         H   il     ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       ⁢     (         y   l     ⁡     (   τ   )       -       y   l     ⁡     (     τ   k     )         )           +       1   2     ⁢       ∑         l   1     ,       l   2     =   0           l   1     ,       l   2     ≠   i           n   -   1       ⁢         (         y   l     ⁡     (   τ   )       -       y     l   1       ⁡     (     τ   k     )         )     T     ⁢       H       l   1     ⁢     l   2         ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       ⁢     (         y     l   1       ⁡     (   τ   )       -       y     l   1       ⁡     (     τ   k     )         )           +     O   ⁡     (              y   ⁡     (   τ   )       -     y   ⁡     (       τ   ^     i     )              3     )               
Dropping the third-order term, the objective function can be approximated as:
 
             F   ⁡     (       y   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )       =       F   ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       +         (       G   i     ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       )     T     ⁢     (         y   i     ⁡     (   τ   )       -       y   i     ⁡     (       τ   ^     i     )         )       +       ∑       l   =   0       l   ≠   i         n   -   1       ⁢         (       G   l     ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       )     T     ⁢     (         y   l     ⁡     (   τ   )       -       y   l     ⁡     (       τ   ^     i     )         )         +       1   2     ⁢       (         y   i     ⁡     (   τ   )       -       y   i     ⁡     (       τ   ^     i     )         )     T     ⁢       H   ii     ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       ⁢     (         y   i     ⁡     (   τ   )       -       y   i     ⁡     (       τ   ^     i     )         )       +         (         y   i     ⁡     (   τ   )       -       y   i     ⁡     (       τ   ^     i     )         )     T     ⁢       ∑       l   =   0       l   ≠   i         n   -   1       ⁢         H   il     ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       ⁢     (         y   l     ⁡     (   τ   )       -       y   l     ⁡     (     τ   k     )         )           +       1   2     ⁢       ∑         l   1     ,       l   2     =   0           l   1     ,       l   2     ≠   i           n   -   1       ⁢         (         y   l     ⁡     (   τ   )       -       y     l   1       ⁡     (     τ   k     )         )     T     ⁢       H       l   1     ⁢     l   2         ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y   i     ⁡     (       τ   ^     i     )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁡     (       τ   ^     i     )         )       ⁢     (         y     l   1       ⁡     (   τ   )       -       y     l   1       ⁡     (     τ   k     )         )                   
The block-diagonal Hessian submatrix H i,i  around the domain point(y 0 (τ k ), . . . , y 0 ({circumflex over (τ)} i ), . . . , y n−1 (τ k ), u({circumflex over (τ)} i )) can then be transformed to a positive definite matrix using a weak Morse transform, as follows.
 
   First, as shown in  FIG. 48 , the Hessian matrix can be decomposed by spectral decomposition as follows: 
           HP   =     P   ⁢           ⁢   Λ                     Hp   α     =         λ   α     ⁢     p   α     ⁢           ⁢   for   ⁢           ⁢   α     =   0       ,   …   ⁢           ,     N   -   1                 Λ   =     (           λ   0         ⋯       0           ⋮       ⋰       ⋮           0       ⋯         λ     N   -   1             )                 H   =     P   ⁢           ⁢   Λ   ⁢           ⁢     P   T             
The N×N matrix P 4802 comprises eigenvectors for the Hessian matrix as columns, and the N×N diagonal matrix Λ  4804  includes the corresponding eigenvalues of the Hessian matrix.
 
   One approach to transforming the Hessian is to create a new matrix {circumflex over (P)}=[{circumflex over (p)} 0 , {circumflex over (p)} 1 , . . . , {circumflex over (p)} n−1 ] in place of the matrix P, where the column vectors are defined as: 
                 P   ^     α     =         {           P   α             if   ⁢           ⁢     λ   α       &gt;   0             0           if   ⁢           ⁢     λ   α       ≤   0           }     ⁢   for   ⁢           ⁢   α     =   0       ,   …   ⁢           ,     N   -   1           
In addition, a new matrix
 
             Λ   ⋒     =     (             λ   ^     0         ⋯       0           ⋮       ⋰       ⋮           0       ⋯           λ   ^       N   -   1             )           
is defined with diagonal entries:
 
                 λ   ^     α     =         {           λ   α             if   ⁢           ⁢     λ   α       &gt;   0             0           if   ⁢             ⁢             ⁢     λ   α       ≤   0           }     ⁢           ⁢   for   ⁢           ⁢   α     =   0       ,   …   ⁢           ,     N   -   1           
Using the previous results, a modified Hessian matrix that is positive definite can be obtained from the modified P matrix and the modified Λ matrices as follows:
 Ĥ={circumflex over (P)}{circumflex over (Λ)}{circumflex over (P)} T   
Furthermore, a reverse for the modified Hessian matrix Ĥ R  can be computed as follows:
 Ĥ R ={circumflex over (P)}{circumflex over (Λ)} R {circumflex over (P)} T   
A reverse is needed because the transformed Hessian may be singular.
 
   A second approach to transforming the Hessian is known as the Strong-Morse Transform, and is used in the Morse-mode forward integration step. In this method, a new modified Λ matrix is produced as follows:
 
{circumflex over (λ)} α =|λ α| for α= 0, . . . , N−1
 
Then, a modified Hessian is computed as follows:
 
Ĥ={circumflex over (P)}{circumflex over (Λ)}{circumflex over (P)} T 
 
This transformed Hessian is non-singular and thus the inverse of the transformed Hessian is well defined, as follows:
 
               H   ^       -   1       =     P   ⁢       Λ   ^       -   1       ⁢     P   T                       Λ   ^     αα     -   1       =     1       λ   ^     α             
In alternative embodiments, the weak Morse transform can be used for points that are outside of the Morse mode, points that are far from a critical point.
 
Using the transformed Hessian, Ĥ i,i , the approximation of the objective function becomes:
 
                       F   ⋓     ⁡     (       y   ⁢           ⁢     (   τ   )       ,     u   ⁢           ⁢     (   τ   )         )       =       ⁢       F   ⁢           ⁢     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁢           ⁢     (       τ   ^               ⁢   i       )         )       +       ⁢                           ⁢       (       G             ⁢   i       ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁢           ⁢     (       τ   ^               ⁢   i       )         )       )     T                     ⁢       (         y             ⁢   i       ⁡     (   τ   )       -       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )         )     +                     ⁢                 ⁢     ∑     l   ⁢           =           ⁢   0               ⁢     n   ⁢           -           ⁢   1                     ⁢     l   ⁢           ≠           ⁢   i         ⁢           ⁢     (       G             ⁢   l       (         y             ⁢   0       ⁢     (     τ             ⁢   k       )     ⁢           ⁢   …     ⁢           ,       y             ⁢   i       ⁢     (               ⁢     τ   ^                 ⁢   i       )       ,     …   ⁢           ⁢       y             ⁢     n   -   1         ⁡     (     τ   k     )         ,   u     ⁢                                         ⁢     (               ⁢     τ   ^                 ⁢   i       )     )     )     T     ⁢     (         y             ⁢   l       ⁡     (   τ   )       -       y             ⁢   l       ⁡     (       τ   ^               ⁢   i       )         )       +       1   2     ⁢       (         y             ⁢   i       ⁡     (   τ   )       -       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )         )     T     ⁢     H   ii                       ⁢     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁢           ⁢     (       τ   ^               ⁢   i       )         )                     ⁢       (         y             ⁢   i       ⁢     (   τ   )       -       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )         )     +       (         y             ⁢   i       ⁢     (   τ   )       -       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )         )     T                       ⁢                 ⁢     ∑     l   ⁢           =           ⁢   0               ⁢     n   ⁢           -           ⁢   1                     ⁢     l   ⁢           ≠           ⁢   i         ⁢       H   il     ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁢           ⁢     (       τ   ^               ⁢   i       )         )                         ⁢       (         y             ⁢   l       ⁡     (   τ   )       -       y             ⁢   l       ⁡     (       τ   ^     k     )         )     +     1   2                       ⁢                 ⁢     ∑       l   1     ,       l   2     ⁢           =           ⁢   0                 ⁢     n   ⁢           -           ⁢   1                     ⁢       l   1     ,       l   2     ⁢           ≠           ⁢   i           ⁢           ⁢       (         y             ⁢     l   1         ⁡     (   τ   )       -       y     l   1       ⁡     (     τ   k     )         )     T     ⁢     H       l   1     ⁢     l   2           ⁢                           ⁢     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁢           ⁢     (       τ   ^               ⁢   i       )         )                     ⁢     (         y     l   2       ⁡     (   τ   )       -       y     l   2       ⁡     (     τ   k     )         )                 
Taking the derivative of this approximate objective function with respect to the i th  partition of the state vector y i (τ) yields the i th  partition of the vector G, as follows:
 
                         ⁢           G             ⋓       i     ⁡     (       y   ⁢           ⁢     (   τ   )       ,     u   ⁢           ⁢     (   τ   )         )       =       ⁢         ∂           ⁢     F             ⋓         ⁢           ⁢     (       y   ⁢           ⁢     (   τ   )       ,           ⁢     u   ⁢           ⁢     (   τ   )         )         ∂           ⁢     y             ⁢   i                         =       ⁢         G   i     ⁡     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁢           ⁢     (       τ   ^               ⁢   i       )         )       +                     ⁢         H   ^     ii     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,   u                           ⁢     (       τ   ^               ⁢   i       )     )     ⁢     (         y             ⁢   i       ⁡     (   τ   )       -       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )         )       +         ∑     l   =   0       n   -   1         l   ≠   i       ⁢           ⁢       H   il     (           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,                           ⁢         y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,     …   ⁢           ⁢       y     n   -   1       ⁡     (     τ   k     )         ,     u   ⁢           ⁢     (       τ   ^               ⁢   i       )         )     ⁢     (         y   l     ⁡     (   τ   )       -       y             ⁢   l       ⁡     (     τ   k     )         )                 
The i th  partition of the vector G evaluated at τ k  is then:
 
                         ⁢           G             ⋓       i     ⁡     (       y   ⁢           ⁢     (   τ   )       ,     u   ⁢           ⁢     (   τ   )         )       =       ⁢       (         ∂           ⁢     F             ⋓         ⁢           ⁢     (       y   ⁢           ⁢     (   τ   )       ,     u   ⁢           ⁢     (   τ   )         )         ∂           ⁢     y             ⁢   i           )     ⁢     |     τ   =     τ   k                         =       ⁢         G   i     (           ⁢           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,           ⁢     …   ⁢           ⁢       y     n   ⁢           -           ⁢   1       ⁡     (     τ   k     )         ,     u   ⁢           ⁢     (       τ   ^               ⁢   i       )         )     +                     ⁢         H   ^     ii     (           ⁢           y   0     ⁡     (     τ   k     )       ⁢           ⁢   …     ⁢           ,       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,           ⁢     …   ⁢           ⁢       y     n   ⁢           -           ⁢   1       ⁡     (     τ   k     )         ,   u                         ⁢       (       τ   ^               ⁢   i       )     )       ⁢           ⁢     (         y             ⁢   i       ⁡     (     τ   k     )       ⁢           -           ⁢       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )         )                 
The Hessian submatrices are therefore given by:
 
                 H             ⋓       ij     ⁡     (       y   ⁢           ⁢     (     τ   k     )       ,     u   ⁢           ⁢     (     τ   k     )         )       =     {                   H   ^     ii     ⁡     (           y             ⁢   0       ⁡     (     τ             ⁢   k       )       ⁢           ⁢   …     ⁢           ,           ⁢       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,           ⁢     …   ⁢           ⁢       y             ⁢     n   ⁢           -           ⁢   1         ⁡     (     τ             ⁢   k       )         ,           ⁢     u   ⁢           ⁢     (       τ   ^               ⁢   i       )         )       ⁢           ⁢   if   ⁢           ⁢   i     ⁢           =           ⁢   j                     H             ⁢   ij       ⁡     (           y             ⁢   0       ⁡     (     τ             ⁢   k       )       ⁢           ⁢   …     ⁢           ,           ⁢       y             ⁢   i       ⁡     (       τ   ^               ⁢   i       )       ,           ⁢     …   ⁢           ⁢       y             ⁢     n   ⁢           -           ⁢   1         ⁡     (     τ             ⁢   _       )         ,           ⁢     u   ⁢           ⁢     (       τ   ^               ⁢   i       )         )       ⁢           ⁢   otherwise     ⁢                           
Again, the off-diagonal Hessian submatrices are not modified by the Strong-Morse Transform since the Hessian is diagonally dominant, and the off-diagonal submatrices contribute very little to the approximation of the objective function. Using assumption 2, above, the block-diagonal Hessian submatrix at τ k  is approximately equal to the block-diagonal Hessian submatrix at {circumflex over (τ)} i     Ĥ   ii ( y   0 (τ k ) . . . ,  y   i ({circumflex over (τ)} i ), . . .  y   n−1 (τ k ),  u ({circumflex over (τ)} i ))≈Ĥ ii ( y   0 ({circumflex over (τ)} i ) . . . ,  y   i ({circumflex over (τ)} i ), . . .  y   n−1 ({circumflex over (τ)} i ),  u ({circumflex over (τ)} i ))=Ĥ ii ( y ({circumflex over (τ)} i ),  u ({circumflex over (τ)} i )) 
Similarly the third assumption, above, implies that
   G   i ( y   0 (τ k ) . . . ,  y   i ({circumflex over (τ)} i ), . . .  y   n−1 (τ k ),  u ({circumflex over (τ)} i ))≈ G   i ( y   0 (τ k ) . . . ,  y   i (τ k ), . . .  y   n−1 (τ k ),  u (τ k ))= G   i ( y (τ k ),  u (τ k ))   H   i ( y   0 (τ k ) . . . ,  y   i ({circumflex over (τ)} i ), . . .  y   n−1 (τ k ),  u ({circumflex over (τ)} i ))≈ H   ij ( y   0 (τ k ) . . . ,  y   i (τ k ), . . .  y   n−1 (τ k ),  u (τ k ))= H   ij ( y (τ k ),  u (τ k )) 
Therefore, the vector G of the approximated objective function {hacek over (F)}(y(τ), u(τ)) evaluated at τ k  is well approximated as follows:
   {hacek over (G)} ( y (τ k ),  u (τ k ))≈ G   i ( y (τ k ),  u (τ k ))+ Ĥ   ii ( y ({circumflex over (τ)} i ),  u ({circumflex over (τ)} i ))( y   i (τ k )− y   i ({circumflex over (τ)} i )) 
and the Hessian submatrix for the approximated objective function evaluated at τ k  can be approximated as
 
                 H             ⋓       ij     ⁡     (       y   ⁢           ⁢     (     τ   k     )       ,     u   ⁢           ⁢     (     τ   k     )         )       =     {                   H   ^     ii     ⁡     (       y   ⁢           ⁢     (       τ   ^               ⁢   i       )       ,     u   ⁢           ⁢     (       τ   ^               ⁢   i       )         )       ⁢           ⁢   if   ⁢           ⁢   i     =   j                     H             ⁢   ij       ⁡     (       y   ⁢           ⁢     (     τ             ⁢   k       )       ,     u   ⁢           ⁢     (     τ   k     )         )       ⁢           ⁢   otherwise     ⁢                           
These approximations then allow for approximate calculation of the w-vector partition and the local ρ for agent i, w i  and τ i , as follows:
 
   
     
       
         
           
             
               
                 
                   
                     w 
                     
                       
                           
                       
                       ⁢ 
                       i 
                     
                   
                   ⁡ 
                   
                     ( 
                     
                       
                         y 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           ( 
                           
                             τ 
                             
                               
                                   
                               
                               ⁢ 
                               k 
                             
                           
                           ) 
                         
                       
                       , 
                       
                         u 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           ( 
                           
                             τ 
                             
                               
                                   
                               
                               ⁢ 
                               k 
                             
                           
                           ) 
                         
                       
                     
                     ) 
                   
                 
                 = 
                   
                 ⁢ 
                 
                   
                     ( 
                     
                       
                         
                           H 
                           ^ 
                         
                         ii 
                       
                       ⁡ 
                       
                         ( 
                         
                           
                             y 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               ( 
                               
                                 
                                   τ 
                                   ^ 
                                 
                                 
                                   
                                       
                                   
                                   ⁢ 
                                   i 
                                 
                               
                               ) 
                             
                           
                           , 
                           
                             u 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               ( 
                               
                                 
                                   τ 
                                   ^ 
                                 
                                 
                                   
                                       
                                   
                                   ⁢ 
                                   i 
                                 
                               
                               ) 
                             
                           
                         
                         ) 
                       
                     
                     ) 
                   
                   
                     - 
                     1 
                   
                 
               
             
           
           
             
               
                   
                 ⁢ 
                 
                   ( 
                   
                     
                       
                         
                           
                               
                           
                           ⁢ 
                           
                             
                               
                                 
                                   G 
                                   
                                     
                                         
                                     
                                     ⋓ 
                                   
                                 
                                 i 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     y 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       ( 
                                       
                                         τ 
                                         
                                           
                                               
                                           
                                           ⁢ 
                                           k 
                                         
                                       
                                       ) 
                                     
                                   
                                   , 
                                   
                                     u 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       ( 
                                       
                                         τ 
                                         
                                           
                                               
                                           
                                           ⁢ 
                                           k 
                                         
                                       
                                       ) 
                                     
                                   
                                 
                                 ) 
                               
                             
                             - 
                           
                         
                       
                     
                     
                       
                         
                           
                             
                               ∑ 
                               
                                 j 
                                 ⁢ 
                                 
                                     
                                 
                                 = 
                                 
                                     
                                 
                                 ⁢ 
                                 0 
                               
                               
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   n 
                                   ⁢ 
                                   
                                       
                                   
                                   - 
                                   
                                       
                                   
                                   ⁢ 
                                   1 
                                 
                               
                             
                             
                               
                                   
                               
                               ⁢ 
                               
                                 j 
                                 ⁢ 
                                 
                                     
                                 
                                 ≠ 
                                 
                                     
                                 
                                 ⁢ 
                                 i 
                               
                             
                           
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             
                               
                                 H 
                                 
                                   
                                       
                                   
                                   ⁢ 
                                   ij 
                                 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     y 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       ( 
                                       
                                         τ 
                                         
                                           
                                               
                                           
                                           ⁢ 
                                           k 
                                         
                                       
                                       ) 
                                     
                                   
                                   , 
                                   
                                     u 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       ( 
                                       
                                         τ 
                                         
                                           
                                               
                                           
                                           ⁢ 
                                           k 
                                         
                                       
                                       ) 
                                     
                                   
                                 
                                 ) 
                               
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               
                                 w 
                                 
                                   
                                       
                                   
                                   ⁢ 
                                   j 
                                 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     y 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       ( 
                                       
                                         τ 
                                         
                                           
                                               
                                           
                                           ⁢ 
                                           
                                             k 
                                             ⁢ 
                                             
                                                 
                                             
                                             - 
                                             
                                                 
                                             
                                             ⁢ 
                                             1 
                                           
                                         
                                       
                                       ) 
                                     
                                   
                                   , 
                                   
                                     u 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       ( 
                                       
                                         τ 
                                         
                                           
                                               
                                           
                                           ⁢ 
                                           
                                             k 
                                             ⁢ 
                                             
                                                 
                                             
                                             - 
                                             
                                                 
                                             
                                             ⁢ 
                                             1 
                                           
                                         
                                       
                                       ) 
                                     
                                   
                                 
                                 ) 
                               
                             
                           
                         
                       
                     
                   
                   ) 
                 
               
             
           
         
       
     
     
       
         
           
             
               ρ 
               i 
             
             ⁡ 
             
               ( 
               
                 
                   y 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     ( 
                     
                       τ 
                       
                         
                             
                         
                         ⁢ 
                         k 
                       
                     
                     ) 
                   
                 
                 , 
                 
                   u 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     ( 
                     
                       τ 
                       
                         
                             
                         
                         ⁢ 
                         k 
                       
                     
                     ) 
                   
                 
               
               ) 
             
           
           = 
           
             
               
                 ( 
                 
                   
                     
                       G 
                       
                         
                             
                         
                         ⋓ 
                       
                     
                     i 
                   
                   ⁡ 
                   
                     ( 
                     
                       
                         y 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           ( 
                           
                             τ 
                             
                               
                                   
                               
                               ⁢ 
                               k 
                             
                           
                           ) 
                         
                       
                       , 
                       
                         u 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           ( 
                           
                             τ 
                             
                               
                                   
                               
                               ⁢ 
                               k 
                             
                           
                           ) 
                         
                       
                     
                     ) 
                   
                 
                 ) 
               
               T 
             
             ⁢ 
             
               
                 w 
                 i 
               
               ⁡ 
               
                 ( 
                 
                   
                     y 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       ( 
                       
                         τ 
                         k 
                       
                       ) 
                     
                   
                   , 
                   
                     u 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       ( 
                       
                         τ 
                         k 
                       
                       ) 
                     
                   
                 
                 ) 
               
             
           
         
       
     
     
       
         
           
             ρ 
             ⁢ 
             
                 
             
             ⁢ 
             
               ( 
               
                 
                   y 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     ( 
                     
                       τ 
                       
                         
                             
                         
                         ⁢ 
                         k 
                       
                     
                     ) 
                   
                 
                 , 
                 
                   u 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     ( 
                     
                       τ 
                       
                         
                             
                         
                         ⁢ 
                         k 
                       
                     
                     ) 
                   
                 
               
               ) 
             
           
           = 
           
             
               ∑ 
               
                 i 
                 = 
                 0 
               
               
                 n 
                 - 
                 1 
               
             
             ⁢ 
             
                 
             
             ⁢ 
             
               
                 ρ 
                 i 
               
               ⁡ 
               
                 ( 
                 
                   
                     y 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       ( 
                       
                         τ 
                         k 
                       
                       ) 
                     
                   
                   , 
                   
                     u 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       ( 
                       
                         τ 
                         k 
                       
                       ) 
                     
                   
                 
                 ) 
               
             
           
         
       
     
     
       
         
           
             
               
                 ⅆ 
                 
                   y 
                   
                     
                         
                     
                     ⁢ 
                     i 
                   
                 
               
               
                 ⅆ 
                 τ 
               
             
             ⁢ 
             
               ( 
               
                 τ 
                 k 
               
               ) 
             
           
           = 
           
             - 
             
               
                 
                   w 
                   i 
                 
                 ⁡ 
                 
                   ( 
                   
                     
                       y 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         ( 
                         
                           τ 
                           k 
                         
                         ) 
                       
                     
                     , 
                     
                       u 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         ( 
                         
                           τ 
                           k 
                         
                         ) 
                       
                     
                   
                   ) 
                 
               
               
                 ρ 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   ( 
                   
                     
                       y 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         ( 
                         
                           τ 
                           k 
                         
                         ) 
                       
                     
                     , 
                     
                       u 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         ( 
                         
                           τ 
                           k 
                         
                         ) 
                       
                     
                   
                   ) 
                 
               
             
           
         
       
     
   
   In step  4630  of  FIG. 46 , the binding variable coefficients may be adjusted. Although the binding variable coefficients might be continuously adjusted, as implied by  FIG. 19 , in at least one embodiment of the present invention, the adjustment is binary and occurs when 
               y   N     -     (       C   ⁢           ⁢     (   y   )       +       ∑     α   =   0         ℳ   2     -   1       ⁢           ⁢       r   α     ⁢     s   α           )       &lt;       u   M     +     ipm_R   ⁢   _switch   ⁢   _offset             
where y N  is the envelope variable z.
 
At the point when the binding variable coefficients are adjusted, they are adjusted to some larger constant value from their original smaller, constant value. At this point, the state vector needs also to be adjusted, so that adjustment of the binding variable coefficient does not force the state vector to an infeasible point. The state vector is adjusted as follows:
 
             y   N     =       (       C   ⁢           ⁢     (   y   )       +       ∑     α   =   0         ℳ   2     -   1       ⁢           ⁢       r   α     ⁢     s   α           )     &lt;       u   M     +     ipm_C   ⁢   _offset2               
where the parameter “ipm_c_offset2” has a value close to one.
 
   The u_and_r module needs to compute a next value for the control vector u during each iteration of the inner loop. This computation is very similar to computation of y vector in y-module agents. As discussed above, the iterative equation for calculation of a next control vector u is:
 
 u   k+1   =u   k −φ k 
 
   For the formulation of the static problem we use a convex quadratic approximation of the Hamiltonian for each iteration. Since we can find an explicit expression for the minimum we do not necessarily need the parameter σ. The u computation only takes one step. 
   Let u k =u| τ     k    be a point of the function H(y,u,p). A second order approximation around u k , given a fixed y k  and p k , using a second order Taylor expansion can be computed as follows: 
                     ⁢     ℋ   ⁢     ❘       y   k     ,     u   k     ,     p   k         ⁢       ≈   ℋ     ⁢     ❘       y   k     ,   u   ,     p   k         +                       ⁢         Ψ   ⁢       |         y   k     ,     u   k     ,     p   k         ⁢     (     u   -     u   k       )       +       1   2     ⁢       (     u   -     u   k       )     T     ⁢       ∂   Ψ       ∂   u           ⁢     ❘       y   k     ,     u   k     ,     p   k         ⁢     (     u   -     u   k       )                   
where Ψ is defined as
 
               ∂   ℋ       ∂   u       .         
This is a quadratic function, but it is not necessarily convex since
 
             ∂   Ψ       ∂   u           
is not necessarily positive definite. A Morse-transform can be used to transform
 
             ∂   Ψ       ∂   u           
to a positive definite convex approximation
 
               ∂     Ψ   ^         ∂   u       ,         
as was done above for the Hessian matrix.
 
   The Hamiltonian matrix can then be approximated as: 
                     ⁢       ℋ   ^     ⁢     ❘       y   k     ,     u   k     ,     p   k         ⁢       ≈   ℋ     ⁢     ❘       y   k     ,     u   k     ,     p   k         +                       ⁢         Ψ   ⁢       |         y   k     ,     u   k     ,     p   k         ⁢     (     u   -     u   k       )       +       1   2     ⁢       (     u   -     u   k       )     T     ⁢       ∂     Ψ   ^         ∂   u           ⁢     ❘       y   k     ,     u   k     ,     p   k         ⁢     (     u   -     u   k       )                             Defining   ⁢           ⁢   Δ   ⁢           ⁢   u     =     u   -     u   k         ,     
     ⁢               ⁢       ℋ   ^     ⁢     ❘       y   k     ,     u   k     ,     p   k         ⁢       ≈   ℋ     ⁢     ❘       y   k     ,     u   k     ,     p   k         +                       ⁢         Ψ   ⁢       |         y   k     ,     u   k     ,     p   k         ⁢   Δ   ⁢           ⁢   u     +       1   2     ⁢           ⁢   Δ   ⁢           ⁢     u   T     ⁢       ∂     Ψ   ^         ∂   u           ⁢     ❘       y   k     ,     u   k     ,     p   k         ⁢     u   ⁢           ⁢   Δ                     
In the following discussion, the dependence on y k , u k , and p k  is not explicitly shown. The minimum of the above function can be found by setting the gradient of the objective function equal to 0 and solving for Δu. Since
 
               ∂     H   ^           ∂   Δ     ⁢           ⁢   u       =       Ψ   +     Δ   ⁢           ⁢     u   T     ⁢       ∂     Ψ   ^         ∂   u           ⁢     |       y   k     ,     u   k     ,     p   k         ⁢     
     ⁢   then                   Δ   ⁢           ⁢   u     =       -       (       ∂     Ψ   ^         ∂   u       )     R       ⁢   Ψ           
where the superscript R denotes the reverse, described above. Defining:
 
             φ   ^     =         (       ∂     Ψ   ^         ∂   u       )     R     ⁢     (   Ψ   )                 then               Δ   ⁢           ⁢   u     =     -     φ   ^             
This allows for iterative update of the control:
   u   k+1   =u   k   +Δu   
However, solving for Ψ and
 
             ∂   Ψ       ∂   u           
is a computationally involved process. The Hamiltonian is expressed as follows:
 
                 ℋ   =     1   -       p   T     ⁢     w   ρ                     =     1   -         ∑     ζ   =   0       N   -   1       ⁢       p   ζ     ⁢     w   ζ         ρ                   
Using the chain rule and the fact that p is treated as a constant during the control computation:
 
             Ψ   α     =         ∂   ℋ       ∂     u   α         =     -       ∑     ζ   =   0       N   -   1       ⁢       p   ζ     ⁡     (         1   ρ     ⁢       ∂     w   ζ         ∂     u   α           -       1     ρ   2       ⁢     w   ζ     ⁢       ∂   ρ       ∂     u   α             )                   
for α=0, . . . , M−1 and
 
               ∂     Ψ   α         ∂     u   β         =           ∂   2     ⁢   H         ∂     u   α       ⁢     ∂     u   β           =     -       ∑     ζ   =   0       N   -   1       ⁢       p   ζ     ⁢           [         1   ρ     ⁢         ∂   2     ⁢     w   ζ           ∂     u   α       ⁢     ∂     u   β             -       1     ρ   2       ⁢       ∂     w   ζ         ∂     u   α         ⁢       ∂   ρ       ∂     u   β           -       1     ρ   2       ⁢     w   ζ     ⁢         ∂   2     ⁢   ρ         ∂     u   α       ⁢     ∂     u   β             +       2     ρ   3       ⁢     w   ζ     ⁢       ∂   ρ       ∂     u   α         ⁢       ∂   ρ       ∂     u   β             ]                     
for α,β=0, . . . , M−1
 
By definition:
 
             w   α     =       ∑     ζ   =   0       N   -   1       ⁢       H     α   ⁢           ⁢   ζ       -   1       ⁢     G   ζ               
for α=0, . . . , N−1, the first derivative of which, with respect to u, is:
 
               ∂     w   α         ∂     u   β         =       ∑     ζ   =   0       N   -   1       ⁢     (           ∂     H   αζ     -   1           ∂     u   β         ⁢     G   ζ       +       H     α   ⁢           ⁢   ζ       -   1       ⁢       ∂     G   ζ         ∂     u   β             )             
for α=0, . . . , N−1 and β=0, . . . , M−1.
 
The inverse of the Hessian matrix H −1  is defined as:
 H −1   H=I   
and can be written in terms of its elements as:
 
               ∑     ζ   =   0       N   -   1       ⁢       H     α   ⁢           ⁢   ζ       -   1       ⁢     H     ζ   ⁢           ⁢   η           =     δ   αη           
for α,β=0, . . . , N−1 and where:
 
             δ     α   ⁢           ⁢   η       =     {         1           if   ⁢           ⁢   α     =   η             0         if   ⁢           ⁢   otherwise                   
Taking the derivative with respect to u on both sides, using the product rule:
 
               ∑     ζ   =   0       N   -   1       ⁢     (           ∂     H   αζ     -   1           ∂     u   γ         ⁢     H   ζη       +       H     α   ⁢           ⁢   ζ       -   1       ⁢         ∂   H     ⁢           ⁢   η       ∂     u   γ             )       =   0         
for α,η=0, . . . , N−1 and γ=0, . . . , M−1. Rearranging terms:
 
               ∑     ζ   =   0       N   -   1       ⁢         ∂     H     α   ⁢           ⁢   ζ       -   1           ∂     u   γ         ⁢     H   ζη         =     -       ∑     η   =   0       N   -   1       ⁢       H     α   ⁢           ⁢   ζ       -   1       ⁢       ∂     H     ζ   ⁢           ⁢   η           ∂     u   γ                     
and multiplying both sides by H −1 :
 
               ∑     ζ   =   0       N   -   1       ⁢       ∑     β   =   0       N   -   1       ⁢         ∂     H     α   ⁢           ⁢   ζ       -   1           ∂     u   γ         ⁢     H     ζ   ⁢           ⁢   η       ⁢     H     η   ⁢           ⁢   β       -   1             =     -       ∑     ζ   =   0       N   -   1       ⁢       ∑     β   =   0       N   -   1       ⁢       H     α   ⁢           ⁢   ζ       -   1       ⁢       ∂     H     ζ   ⁢           ⁢   η           ∂     u   γ         ⁢     H     η   ⁢           ⁢   β       -   1                     
the expression can be simplified to:
 
               ∑     ζ   =   0       N   -   1       ⁢         ∂     H     α   ⁢           ⁢   ζ       -   1           ∂     u   γ         ⁢     δ     ζ   ⁢           ⁢   β           =       ∂     H     α   ⁢           ⁢   β       -   1           ∂     u   γ               
providing an explicit expression for the elements of H −1   
               ∂     H     -   1           ∂     u   r         =     -       ∑     η   =   0       N   -   1       ⁢       ∑     β   =   0       N   -   1       ⁢       H   αζ     -   1       ⁢       ∂     H   ζη         ∂     u   r         ⁢     H   ηβ     -   1                     
Substituting the expression for H −1  into:
 
               ∂     w   α         ∂     u   β         =       ∑     ζ   =   0       N   -   1       ⁢     (           ∂     H   αζ     -   1           ∂     u   β         ⁢     G   ζ       +       H   αζ     -   1       ⁢       ∂     G   ζ         ∂     u   β             )             
provides:
 
                     ∂     w   α         ∂     u   β         =       ⁢       ∑     ζ   =   0       N   -   1       ⁢     (           ∂     H   αζ     -   1           ∂     u   β         ⁢     G   ζ       +       H   αζ     -   1       ⁢       ∂     G   ζ         ∂     u   β             )                   =       ⁢       ∑     ζ   =   0       N   -   1       ⁢     (         (     -       ∑     η   =   0       N   -   1       ⁢       ∑     ω   =   0       N   -   1       ⁢       H   αη     -   1       ⁢       ∂     H   ηω         ∂     u   β         ⁢     H   ωζ     -   1               )     ⁢     G   ζ       +       H   αζ     -   1       ⁢       ∂     G   ζ         ∂     u   β             )                   =       ⁢       -       ∑     ζ   =   0       N   -   1       ⁢       ∑     η   =   0       N   -   1       ⁢       ∑     ω   =   0       N   -   1       ⁢       H   αη     -   1       ⁢       ∂     H   ηω         ∂     u   β         ⁢     H   ωζ     -   1       ⁢     G   ζ               +       ∑     η   =   0       N   -   1       ⁢       H   αη     -   1       ⁢       ∂     G   η         ∂     u   β                           =       ⁢       -       ∑     η   =   0       N   -   1       ⁢       ∑     ω   =   0       N   -   1       ⁢       H   αη     -   1       ⁢       ∂     H   ηω         ∂     u   β         ⁢     w   ω             +       ∑     η   =   0       N   -   1       ⁢       H   αη     -   1       ⁢       ∂     G   η         ∂     u   β                           =       ⁢       ∑     η   =   0       N   -   1       ⁢       H   αη     -   1       ⁡     (         ∂     G   η         ∂     u   β         -       ∑     ω   =   0       N   -   1       ⁢         ∂     H   ηω         ∂     u   β         ⁢     w   ω           )                     
for α=0, . . . , N−1 and β=0, . . . , M−1.
 
The second derivative is
 
                 ∂   2     ⁢     w   α           ∂     u   β       ⁢     ∂     u   γ           =       ∑     η   =   0       N   -   1       ⁢     {           ∂     H   αη     -   1           ∂     u   γ         ⁢     (         ∂     G   η         ∂     u   β         -       ∑     ω   =   0       N   -   1       ⁢         ∂     H   ηω         ∂     u   β         ⁢     w   ω           )       -       H   αη     -   1       ⁡     (       ∑     ω   =   0       N   -   1       ⁢         ∂     H   ηω         ∂     u   β         ⁢       ∂     w   ω         ∂     u   γ             )         }             
Note that the expression is simplified by the fact that the second derivative of the gradient of F and Hessian matrix with respect to u is 0 since these terms are linear in u. Substituting the expression for H −1  into the above expression provides:
 
                 ∂   2     ⁢     w   α           ∂     u   β       ⁢     ∂     u   γ           =       -       ∑     ζ   =   0       N   -   1       ⁢     {       H   αζ     -   1       ⁢       ∑     θ   =   0       N   -   1       ⁢         ∂     H   ζθ         ∂     u   γ         ⁢       ∑     η   =   0       N   -   1       ⁢       H   θη     -   1       ⁡     (         ∂     G   η         ∂     u   β         -       ∑     ω   =   0       N   -   1       ⁢         ∂     H   ηω         ∂     u   β         ⁢     w   ω           )               }         -       ∑     η   =   0       N   -   1       ⁢       H   αη     -   1       ⁡     (       ∑     ω   =   0       N   -   1       ⁢         ∂     H   ηω         ∂     u   β         ⁢       ∂     w   ω         ∂     u   γ             )                 
This can be simplified as follows:
 
                             ⁢       ∂   2     ⁢           ⁢     w             ⁢   α               ∂           ⁢     u             ⁢   β         ⁢           ⁢     ∂           ⁢     u             ⁢   γ             =       ⁢       -       ∑     ζ   ⁢           =           ⁢   0               ⁢     N   ⁢           -           ⁢   1         ⁢       H             ⁢   αζ       -   1       ⁢     (       ∑     θ   ⁢           =           ⁢   0               ⁢     N   ⁢           -           ⁢   1         ⁢         ∂           ⁢     H             ⁢   ζθ           ∂           ⁢     u             ⁢   γ           ⁢       ∂           ⁢     w             ⁢   θ           ∂           ⁢     u             ⁢   β               )           -                     ⁢       ∑     η   ⁢           =           ⁢   0               ⁢     N   ⁢           -           ⁢   1         ⁢       H             ⁢   αη       -   1       ⁡     (       ∑     ω   ⁢           =           ⁢   0               ⁢     N   ⁢           -           ⁢   1         ⁢         ∂           ⁢     H             ⁢   ηω           ∂           ⁢     u             ⁢   β           ⁢       ∂           ⁢     w             ⁢   ω           ∂           ⁢     u             ⁢   γ               )                     =       ⁢       -       ∑     ζ   =   0       N   -   1       ⁢       H   αζ     -   1       ⁡     (       ∑     θ   =   0       N   -   1       ⁢         ∂     H   ζθ         ∂     u   γ         ⁢       ∂     w   θ         ∂     u   β             )           -       ∑     ζ   =   0       N   -   1       ⁢       H   αζ     -   1       ⁡     (       ∑     θ   =   0       N   -   1       ⁢         ∂     H   ζθ         ∂     u   β         ⁢       ∂     w   θ         ∂     u   γ             )                       =       ⁢     -       ∑     ζ   =   0       N   -   1       ⁢       H   αζ     -   1       ⁡     (         ∑     θ   =   0       N   -   1       ⁢         ∂     H   ζθ         ∂     u   γ         ⁢       ∂     w   θ         ∂     u   β             +         ∂     H   ζθ         ∂     u   β         ⁢       ∂     w   θ         ∂     u   γ             )                       
Next the denominator, ρ, is defined as
 
           ρ   =       ∑     ζ   =   0       N   -   1       ⁢       G   ζ     ⁢     w   ζ               
The first and second derivatives with respect to u are, using the product rule,
 
               ∂   ρ       ∂     u   α         =       ∑     ζ   =   0       N   -   1       ⁢     (           ∂     G   ζ         ∂     u   α         ⁢     w   ζ       +       G   ζ     ⁢       ∂     w   ζ         ∂     u   α             )             
for α=0, . . . , M−1
 
and
 
                 ∂   2     ⁢   ρ         ∂     u   α       ⁢     ∂     u   β           =       ∑     ζ   =   0       N   -   1       ⁢     (           ∂     G   ζ         ∂     u   α         ⁢       ∂     w   ζ         ∂     u   β           +         ∂     G   ζ         ∂     u   β         ⁢       ∂     w   ζ         ∂     u   α           +       G   ζ     ⁢         ∂   2     ⁢     w   ζ           ∂     u   α       ⁢     u   β             )             
for α, β=0, . . . , M−1
 
Next, a formal distributed algorithm is discussed. Multiplying both sides of:
 
             φ   ^     =         (       ∂     Ψ   ^         ∂   u       )     R     ⁢     (   Ψ   )                   by   ⁢     :                   (       (       ∂     Ψ   ^         ∂   u       )     R     )       -   1           
results in:
 
                 ∂     Ψ   ^         ∂   u       ⁢     φ   ^       =   Ψ         
Expressing the relationship in terms of the r th  agent:
 
                 ∂       Ψ   ^     r         ∂     u   r         ⁢       φ   ^     r       +       ∑       s   =   0       s   ≠   r         M   -   1       ⁢         ∂     Ψ   rs         ∂     u   s         ⁢       φ   ^     s     ⁢     Ψ   r               
for i=0, . . . , M−1. Rearranging terms:
 
                 ∂       Ψ   ^     r         ∂     u   r         ⁢       φ   ^     r       =       Ψ   r     -       ∑       s   =   0       s   ≠   r         M   -   1       ⁢         ∂     Ψ   rs         ∂     u   s         ⁢       φ   ^     s                 
Using the notation:
 
               Ψ   ~     r     =     (       Ψ   r     -       ∑       s   =   0       s   ≠   r         M   -   1       ⁢         ∂     Ψ   rs         ∂     u   s         ⁢       φ   ^     s           )           
Multiplying both sides by the inverse of
 
             ∂       Ψ   ^     r         ∂     u   r             
provides:
 
               φ   ^     r     =         (       ∂       Ψ   ^     r         ∂     u   r         )     R     ⁢       Ψ   ~     r             
results in:
 Δ u   r =−{circumflex over (φ)} r   
Thus, the control can be updated by:
   u   k+1   r   +=u   k   r   +Δu   r   
   The definition for the Hamiltonian is now written in terms of agents: 
           ℋ   =     1   -         ∑     i   =   0       n   -   1       ⁢       ∑     ζ   =   0         N   ′     -   1       ⁢       p   ζ   i     ⁢     w   ζ   i           ρ                 Then               Ψ   α   r     =         ∂   ℋ       ∂     u   α   r         =     -       ∑     i   =   0       n   -   1       ⁢       ∑     ζ   =   0         N   ′     -   1       ⁢       p   ζ   i     ⁡     [         1   ρ     ⁢       ∂     w   ζ   i         ∂     u   α   r           -       1     ρ   2       ⁢     w   ζ   i     ⁢       ∂   ρ       ∂     u   α   r             ]                         and                 ∂     Ψ   α   r         ∂     u   β   s         =     -       ∑     i   =   0       n   -   1       ⁢       ∑     ζ   =   0         N   i     -   1       ⁢       p   ζ   i     ⁢           ⁢     
     [         1   ρ     ⁢         ∂   2     ⁢     w   ζ   i           ∂     u   α   r       ⁢     ∂     u   β   s             -       1     ρ   2       ⁢       w   ζ   i       ∂     u   β   s         ⁢       ∂   ρ       ∂     u   α   r           -       1     ρ   2       ⁢       ∂     w   ζ   i         ∂     u   α   s         ⁢       ∂   ρ       ∂     u   β   r           -       1     ρ   2       ⁢     w   ζ   i     ⁢         ∂   2     ⁢   ρ         ∂     u   α   r       ⁢     ∂     u   β   s             +       2     ρ   3       ⁢     w   ζ   i     ⁢       ∂   ρ       ∂     u   α   r         ⁢       ∂   ρ       ∂     u   β   s             ]                 
Considering that:
 
             w   α   i     =       ∑     ζ   =   0         N   ′     -   1       ⁢         (     H   ii     )     αζ     -   1       ⁢       G   ~     ζ   i               
for α=0, . . . , n i  and that:
 
               G   ~     ζ   i     =       G   ζ   i     -       ∑       j   =   0       j   ≠   i         n   -   1       ⁢       ∑     η   =   0         N   j     -   1       ⁢         (     H   ij     )       ζ   ⁢           ⁢   η       ⁢     w   η   j                   
each element a of the w vector for agent I can be expressed as:
 
               ∂     w   α   i         ∂     u   β   r         =       ∑     ζ   =   0         N   i     -   1       ⁢     {           ∂       (     H   ii     )     αζ     -   1           ∂     u   β   r         ⁢       G   ~     ζ   i       +         (     H   ii     )       α   ⁢           ⁢   ζ       -   1       ⁢       ∂       G   ~     ζ   i         ∂     u   β   r             }             
Substituting the definition of H −1  in the previous expression:
 
               ∂     w   α   i         ∂     u   β   r         =       ∑     ζ   =   0         N   i     -   1       ⁢           ⁢     {           (     -       ∑     η   =   0         N   i     -   1       ⁢       ∑     ω   =   0         N   i     -   1       ⁢         (     H   ii     )     αη     -   1       ⁢       ∂     H   ηω   ii         ∂     u   β         ⁢       (     H   ii     )     ωζ     -   1               )                   G   ~     ζ   i     +         (     H   ii     )     αζ     -   1       ⁢       ∂       G   ~     ζ   i         ∂     u   β   r                   }             
Rearranging terms:
 
               ∂     w   α   i         ∂     u   β   r         =         ∑     ζ   =   0         N   i     -   1       ⁢         (     H   ii     )     αζ     -   1       ⁢       ∂       G   ~     ζ   i         ∂     u   β   r             -       ∑     η   =   0         N   i     -   1       ⁢       ∑     ω   =   0         N   i     -   1       ⁢         (     H   ii     )     αη     -   1       ⁢       ∂     H   ηω   ii         ∂     u   β   r         ⁢     w   ω   i                             ∂     w   α   i         ∂     u   β   r         =       ∑     ζ   =   0         N   i     -   1       ⁢         (     H   ii     )     αζ     -   1       ⁢     {         ∂       G   ~     ζ   i         ∂     u   β   r         -       ∑     ω   =   0         n   i     -   1       ⁢         ∂       (     H   ii     )     ζω         ∂     u   β   r         ⁢     w   ω   i           }               
Since the second derivatives of G and H with respect to u are 0:
 
                 ∂   2     ⁢     w   α   i           ∂     u   β   r       ⁢     ∂     u   γ   s           =         ∑     ζ   =   0         N   i     -   1       ⁢         ∂       (     H   ii     )     αζ     -   1           ∂     u   γ   s         ⁢     {         ∂       G   ~     ζ   i         ∂     u   β   r         -       ∑     ω   =   0         n   i     -   1       ⁢         ∂     H   ζω   ii         ∂     u   β         ⁢     w   ω   i           }         -         (     H   ii     )     αζ     -   1       ⁢     (       ∑     ω   =   0         n   i     -   1       ⁢         ∂     H   ζω   ii         ∂     u   β   r         ⁢       ∂     w   ω   i         ∂     u   γ   s             )       +         (     H   ii     )     αζ     -   1       ⁢     (         ∂   2     ⁢       G   ~     ζ   i           ∂     u   β   r       ⁢     ∂     u   γ   s           )               
However:
 
                 ∂   2     ⁢     w   α   i           ∂     u   β   r       ⁢     ∂     u   γ   s           =         ∑     ζ   =   0         N   i     -   1       ⁢     (     -       ∑     κ   =   0         N   i     -   1       ⁢       ∑     θ   =   0         N   i     -   1       ⁢         (     H   ii     )     ακ     -   1       ⁢       ∂     H   κθ   ii         ∂     u   γ   s         ⁢       ∂     w   θ   i         ∂     u   β   r                 )       -         (     H   ii     )     αζ     -   1       ⁢     (       ∑     ζ   =   0         N   i     -   1       ⁢         ∂     H   ζω   ii         ∂     u   β   r         ⁢       ∂     w   ω   i         ∂     u   γ   s             )       +         (     H   ii     )     αζ     -   1       ⁢     (         ∂   2     ⁢       G   ~     ζ   i           ∂     u   β   r       ⁢     ∂     u   γ   s           )               
Separating terms and changing indices of the last two terms:
 
                       ∂   2     ⁢     w   α   i           ∂     u   β   r       ⁢     ∂     u   γ   s           =       ⁢       (     -       ∑     κ   =   0         N   i     -   1       ⁢       ∑     θ   =   0         N   i     -   1       ⁢         (     H   ii     )     ακ     -   1       ⁢       ∂     H   κθ   ii         ∂     u   γ   s         ⁢       ∂     w   θ   i         ∂     u   β   r                 )     -                     ⁢       (       ∑     ω   =   0         N   i     -   1       ⁢         (     H   ii     )     αζ     -   1       ⁢       ∑     ω   =   0         n   i     -   1       ⁢         ∂     H   ζω   ii         ∂     u   β   r         ⁢       ∂     w   ω   i         ∂     u   γ   s                 )     +       ∑     ζ   =   0         N   i     -   1       ⁢         (     H   ii     )     αζ     -   1       ⁢     (         ∂   2     ⁢       G   ~     ζ   i           ∂     u   β   r       ⁢     ∂     u   γ   s           )                       =       ⁢       (     -       ∑     κ   =   0         N   i     -   1       ⁢       ∑     θ   =   0         N   i     -   1       ⁢         (     H   ii     )     ακ     -   1       ⁢       ∂     H   κθ   ii         ∂     u   γ   s         ⁢       ∂     w   θ   i         ∂     u   β   r                 )     -                     ⁢       (       ∑     κ   =   0         N   i     -   1       ⁢         (     H   ii     )     ακ     -   1       ⁢       ∑     θ   =   0         N   i     -   1       ⁢         ∂     H   κθ   ii         ∂     u   β   r         ⁢       ∂     w   θ   i         ∂     u   γ   s                 )     +                     ⁢       ∑     κ   =   0         N   i     -   1       ⁢         (     H   ii     )     ακ     -   1       ⁢     (         ∂   2     ⁢       G   ~     κ   i           ∂     u   β   r       ⁢     ∂     u   γ   s           )                     =       ⁢     -       ∑     κ   =   0         N   i     -   1       ⁢         (     H   ii     )     ακ     -   1       ⁢     {         ∑     θ   =   0         N   i     -   1       ⁢     (                 ∂     H   κθ   ii         ∂     u   γ   s         ⁢       ∂     w   θ   i         ∂     u   β   r           +                   ∂     H   κθ   ii         ∂     u   β   r         ⁢       ∂     w   θ   i         ∂     u   γ   s                 )       -     (         ∂   2     ⁢       G   ~     κ   i           ∂     u   β   r       ⁢     ∂     u   γ   s           )       }                       
Therefore:
 
               ∂       G   ~     α   i         ∂     u   β   r         =         ∂     G   α   i         ∂     u   β   r         -         ∑     j   =   0       n   -   1         j   ≠   i       ⁢           ⁢       ∑     γ   =   0         N   j     -   1       ⁢           ⁢     (           ∂     H   αζ   ij         ∂     u   β   r         ⁢     w   ζ   j       +       H   αζ   ij     ⁢       ∂     w   ζ   j         ∂     u   β   r             )                 
and using the fact that the second derivatives of G and H with respect to u are 0:
 
                 ∂   2     ⁢       G   ~     α   i           ∂     u   β   r       ⁢     ∂     u   γ   s           =     -         ∑     j   =   0       n   -   1         j   ≠   i       ⁢           ⁢       ∑     γ   =   0         N   j     -   1       ⁢     (           ∂     H   αζ   ij         ∂     u   β   r         ⁢       ∂     w   ζ   j         ∂     u   γ   s           +         ∂     H   αζ   ij         ∂     u   γ   s         ⁢       ∂     w   ζ   j         ∂     u   β   r           +       H   αζ   ij     ⁢         ∂   2     ⁢     w   ζ   j           ∂     u   β   r       ⁢     ∂     u   γ   s               )                 
Finally:
 
           ρ   =       ∑     i   =   0       n   -   1       ⁢           ⁢       ∑     ζ   =     1   -   0           N   j     -   1       ⁢           ⁢       G   ζ   i     ⁢     w   ζ   i                 
and its derivatives are:
 
   
     
       
         
           
             
               ∂ 
               ρ 
             
             
               ∂ 
               
                 u 
                 α 
                 r 
               
             
           
           = 
           
             
               
                 ∑ 
                 
                   i 
                   = 
                   0 
                 
                 
                   n 
                   - 
                   1 
                 
               
               ⁢ 
               
                   
               
               ⁢ 
               
                 
                   ∑ 
                   
                     ζ 
                     = 
                     0 
                   
                   
                     
                       N 
                       i 
                     
                     - 
                     1 
                   
                 
                 ⁢ 
                 
                   
                     
                       ∂ 
                       
                         G 
                         ζ 
                         i 
                       
                     
                     
                       ∂ 
                       
                         u 
                         α 
                         r 
                       
                     
                   
                   ⁢ 
                   
                     w 
                     ζ 
                     i 
                   
                 
               
             
             + 
             
               
                 G 
                 ζ 
                 i 
               
               ⁢ 
               
                 
                   
                     ∂ 
                     2 
                   
                   ⁢ 
                   
                     w 
                     ζ 
                     i 
                   
                 
                 
                   ∂ 
                   
                     u 
                     α 
                     r 
                   
                 
               
             
           
         
       
     
     
       
         
           
             
               
                 ∂ 
                 2 
               
               ⁢ 
               ρ 
             
             
               
                 ∂ 
                 
                   u 
                   α 
                   r 
                 
               
               ⁢ 
               
                 ∂ 
                 
                   u 
                   β 
                   s 
                 
               
             
           
           = 
           
             
               
                 ∑ 
                 
                   i 
                   = 
                   0 
                 
                 
                   N 
                   - 
                   1 
                 
               
               ⁢ 
               
                 
                   ∑ 
                   
                     ζ 
                     = 
                     0 
                   
                   
                     
                       N 
                       i 
                     
                     - 
                     1 
                   
                 
                 ⁢ 
                 
                   
                     
                       ∂ 
                       
                         G 
                         ζ 
                         i 
                       
                     
                     
                       ∂ 
                       
                         u 
                         α 
                         r 
                       
                     
                   
                   ⁢ 
                   
                     
                       ∂ 
                       
                         w 
                         ζ 
                         i 
                       
                     
                     
                       ∂ 
                       
                         u 
                         β 
                         s 
                       
                     
                   
                 
               
             
             + 
             
               
                 
                   ∂ 
                   
                     G 
                     ζ 
                     i 
                   
                 
                 
                   ∂ 
                   
                     u 
                     α 
                     s 
                   
                 
               
               ⁢ 
               
                 
                   ∂ 
                   
                     w 
                     ζ 
                     i 
                   
                 
                 
                   ∂ 
                   
                     u 
                     β 
                     r 
                   
                 
               
             
             + 
             
               
                 G 
                 ζ 
                 i 
               
               ⁢ 
               
                 
                   
                     ∂ 
                     2 
                   
                   ⁢ 
                   
                     w 
                     ζ 
                     i 
                   
                 
                 
                   
                     ∂ 
                     
                       u 
                       α 
                       r 
                     
                   
                   ⁢ 
                   
                     ∂ 
                     β 
                     s 
                   
                 
               
             
           
         
       
     
   
   The propagation of the costate variables is different from that of state and control variables, which are tightly coupled and proceed together forward in time. In contrast, the costate variables are propagated backwards in time, independently of the integration of the state and control variables. In order to do this, a database of values calculated in they module and u_and_r module is needed. Since the differential equation of p is linear in p, its solution can be approximated by an interpolation method rather by integration. For many applications the value of the vector p is almost constant with the exception of the transient around the convergence horizon. Also the matrix A in the equation for p does not depend on p. 
   To compute p, a set of the node values of y and u trajectories calculated by the y module and u_and_r module is created for use by the p module. The terminal conditions p(T) are computed first, and then earlier values of p are calculated in a loop which interpolates backwards to the beginning of the y computation. 
   
     
       
         
           
             
               ⅆ 
               p 
             
             
               ⅆ 
               τ 
             
           
           = 
           AP 
         
       
     
     
       
         
           A 
           = 
           
             
               
                 1 
                 ρ 
               
               ⁢ 
               
                 
                   ( 
                   
                     
                       ∂ 
                       w 
                     
                     
                       ∂ 
                       y 
                     
                   
                   ) 
                 
                 T 
               
             
             - 
             
               
                 1 
                 
                   ρ 
                   2 
                 
               
               ⁢ 
               
                 
                   ∂ 
                   ρ 
                 
                 
                   ∂ 
                   y 
                 
               
               ⁢ 
               
                 w 
                 T 
               
             
           
         
       
     
   
   For practicality, we divide the problem up into increments of size τ k−1 −τ k  where τ k−1 &lt;τ k  where the matrix A is assumed to be nearly constant. Then, integrating backwards the solution at τ k−1  given τ k  is
 
 p (τ k−1 )=exp( A (τ k−1 −τ k )) p (τ k )
 
The exponential of this equation is solved by the approximation
 
             exp   ⁡     (     A   ⁡     (       τ     κ   -   1       -     τ   κ       )       )       ≅       ∑     ζ   =   0       N   -   1       ⁢           A   ζ     ⁡     (       τ     κ   -   1       -     τ   κ       )       ζ       ζ   !               
where N is the degree of the approximation.
 
   Repair 
   So far, the discussion has focused on specific computational methods for generic optimization. However, an additionla aspect of real-world optimization concerns repeated optimization over real periods of time. For example, an enterprise, such as a manufacturing facility, may employ optimization techniques for optimizing all the various aspects of running the manufacturing plant, including sheduling of employees, ordering of parts and raw materials, distrubution of parts, raw materials, and intermediate manufactured articles within the manufacturing plant, and delivery of product to customers, among other things. In general, in the extremely high dimensional problem domains for such enterprise optimizations, and because of the many uncontrolled events and changes that arise during operation of a manufacturing plant, optimization may need to be carried out repeatedly, over time, in order to readjust the control of the manufacturing plant to respond to changing conditions. Similar considerations apply to many other computational tasks, in additio to optimization. for example, modelling of a complex system may require intermediate computation of system state in order to respond to various randomly-generated events and perturbations to the system, introduced so that the model reflects real time response of the system to a changing enviroment. 
     FIG. 49  illustrates this repeated optimization strategy. In  FIG. 49 , and in subsequent FIGS.  50  and  52 – 56 , the horizontal axis  4902  corresponds to time. In  FIG. 49 , a near-optimal state function for a system, including state variables and control variables, is computed via optimization and/or modeling techniques, at each point in time represented in  FIG. 49  by a vertical dashed line, such as vertical dashed line  4904 . The computed state functions are represented in  FIG. 49  as curved segments, such as curved segment  4906 , that begins at a point in time  4908  at which computation of the state function begins ans extends to another point in time  4901 . Of course, computed state fucntions are, in real-world applications, complex control trajectories in hyper-dimensional spaces, but two-dimensional control trajectories, such as curved segment  4906 , are used for clarity of illustration. The diagonal, dashed lines, such as diagonal dash line  4912 , represents the computation of the state function via integration with respect to a continuation iteration variable τ, described above, and takes some amount of time  4914  related to the computation horizon T for the optimization technique by which the state function is calculated. When a new state function is calculated, as, for example, new state function  4916  applicable to the range of time beginning at time  4918  and ending at time  4920 , the new state function likely relatively closely overlaps, but does not coincide with the previous state function, in the present case state function  4906 . Thus, there is a discrepancy  4922  or discontinuity between successively calculated state functions. This discontinuity arises from a change in conditions and environment of the system between the time the previous state function was calculated and the time that the current state function calculation began, in this case the events and changes that occur in time interval beginning at time  4908  and ending at time  4918 . 
   Unfortunately, these discontinuities represent disturbances in the system. For example, in the manufacturing plant, a discontinuity may reflect a change in scheduling of workers, a change in departure times and destinations for delivery trucks, and other such changes. However, complex systems tend to have inertia, and major disturbances are difficult to overcome, may take finite amounts of time and resources to overcome that themselves change the environment of the system to a degree that the previously calculated optimal state function is no longer optimal, and may cause ripples and cascading effects that further perturb a complex system. 
   One possible approach to ameliorating the discontinuities arising from repeated optimization is to shorten the time between optimizations.  FIG. 50  illustrates a shortened interval between optimizations for the state function calculations illustrated in  FIG. 49 . By more efficiently computing the state functions, represented in  FIG. 50  by increased slope of the dashed lines, such as dashed line  5002 , representing state function calculation, the optimizations may be carried out in a shorter amount of time ΔT  5004 . However, the effect of doing so increases the frequency of disturbance, and the increased frequency of disturbance may itself create rippling and cascading effects that eventually lead to greater and greater instabilities. In certain cases, even greater discrepancies, such as discrepancy  5006 , may arise as a result of underdamping of the effects of the occurrence of certain kind of events in very short time intervals. 
     FIG. 51  illustrates an analogy for a second consequence arising from repeated optimizations. Note that, in  FIGS. 49–50 , it is assumed that the state functions are calculated anew without regard for the previously calculated state functions, other than at the point where the previously calculated state function is abandoned and a new state function begins to control the system. However, failure to take into account the history of previous state function calculations may result in state function solutions that are exceedingly non-optimal over longer periods of time. As an analogy, consider the task of traversing terrain represented by a simple topographic map, in  FIG. 51 , from point A  5102  to point B  5104 . A first path  5106  from point A towards point B represents a path that might be followed by an individual relying only on a compass, who frequently readjusts his or her trajectory based only on an instantaneous compass reading without even waiting for the compass needle motion to settle. Such a path represents an exceedingly non-optimal solution, in time and distance, for traversing the terrain between points A and B. 
   A second path  5108  is a path that might be taken by an individual who recognizes that the compass must be allowed to settle, before reading, who recognizes that a compass may be perturbed by local conditions, such as an outcropping of magnetite in a cliff face, and who recognizes that if an obstacle is encountered, it may be prudent to follow that obstacle rather than relying on the compass for direction. Obviously, the second path  5108  represents a far more efficient and effective strategy for traversing the terrain from point A to point B. 
   A third path  5110  represents a path taken by an individual capable of reading a topographic math, properly using a compass, who recognizes that a bridge across a river is a worthwhile intermediate point to aim for, and who recognizes that it is far easier to walk along a path through the forest than to attempt to blaze a new trail. In the analogy of  FIG. 51 , the individual who can rely on past information, rely on information in addition to a simple, limited current state, and who can project future solutions based on additional information is the individual who is most likely to find an optimal path between point A and point B. Returning to the repeated optimization problem illustrated in  FIGS. 49–50 , it is desirable to find a way to employ previously calculated state functions in order to employ more information over longer periods of time, rather than to simply consider local, instantaneous state information. 
   As a result of the considerations described above with reference to  FIGS. 49–51 , an approach to a more or less continuous optimization, over time, using historical information and updated projections, represents an additional feature of embodiments of the present invention directed towards optimization, referred to as a “repair-based continuous optimization approach” (“RCO”). 
     FIG. 52  illustrates the concept of a sliding window. As shown in  FIG. 52 , a prior window  5202  beginning at time t 1 −ΔT  5204  and ending at time t 1 +T−Δτ  5206 , of time-length T, is used as a frame of reference for computing a state function. Since that time, an amount of time ΔT  5208  has elapsed. Therefore, a new, current window  5210  of length T in time provides the frame of reference at the current time t 1    5212 . 
   As shown in  FIG. 53 , the state function calculated in the prior window  5302  extends from time t 1 −ΔT  5204  to time t 1 −ΔT+T  5206 .  FIG. 54  illustrates extension of the previously calculated state function to extend over the current window. As shown in  FIG. 54 , the previously calculated state function  5303  is extended  5402  to the end of the current window at time t 1 +T  5404 .  FIG. 55  illustrates consideration, at the current point of time, of past events and environmental changes that may result in the previously calculated state function being out-of-date, and not optimal, at the current point of time, as well as an approach to incorporating information from the previous state function in finding a new, optimal state function for the current frame of reference represented by the current window. 
   In  FIG. 55 , vertical lines  5502 – 5505  represent events that occurred during the time of calculation of the state function, t 1 −ΔT, and the current time, t 1 . Such events may alter the state and environment of the system in such a way that the previously calculated state function is no longer optimal. For example, in the manufacturing plant, a large number of employees may have decided to go on strike, or a flu virus may have greatly increased the number of workers unavailable for being scheduled to perform tasks within the manufacturing plant. Therefore, previously calculated optimal scheduling may be exceedingly suboptimal in the face of a shortage of workers. As shown in  FIG. 55 , various projections of time-dependent parameters on which the state function was calculated may have changed at time t 1 . The change in these projections is represented by the dashed line  5506 . For example, in the context of the manufacturing plant, a major customer may have indicated an immediate need for twice as many products as the customer normally buys, or the arrival of a new product produced by a competitor may suddenly dampen the demand for the products produced by the manufacturing plant. Thus, previously projected demand for products may not reflect the current and forecast demand. 
   In addition, priorities originally assigned to various considerations and parameters on which the state function is based may have changed. The dotted line  5508  in  FIG. 55  indicates a change in priorities that has arisen for the current frame of reference with respect to a prior frame of reference. For example, at the time that the state function was calculated, there was no knowledge of an impending road closure on Thursday, announced during the time interval t 1 −ΔT to t 1 . Knowledge of the road closure on Thursday may greatly increase the priority of scheduling south-side-resident workers for work at times when the bridge closure will least affect their commutes. 
   Also, in order to prevent serious disturbances to the operation of a system, the system may require that previously calculated state for the system needs to be adhered to, regardless of whether or not a newly calculated state function would provide for the same state over certain intervals. For example, in  FIG. 55 , two time windows  5510  and  5512 , beginning at times T 2  and T 4    5514  and  5516 , respectively, and ending at times T 3  and T 5    5518  and  5520 , respectively, have been established, within which the old value of the state function  5522  and  5524  needs to be maintained in a newly computed state function. Thus, historical information obtained in the previous state-function calculation is maintained in the current frame of reference. Such fence intervals may be necessary to prevent undue disruption and disturbances in system control established based on the previously calculated state function. 
   Using the knowledge of events, knowledge of changes in circumstances and priorities, and the need for maintaining certain state previously calculated, the previously calculated and extended state function can be modified, without a de novo state function calculation, in order to produce a near-optimal state function for the current time window.  FIG. 56  shows repair of the previously calculated state function in order to provide a near optimal state function for a current time window. As shown in  FIG. 56 , the repair of the previously calculated state function involves computing differences, such as difference  5602 , at points in time within the current time window, between the previously calculated state function and a new, near optimal, repaired state function for the current time window. Computation of the differences, such as difference  5602 , takes into account the occurrence of events in the prior interval of time t 1 −ΔT  5208 , changes in circumstances and priorities ( 5508  and  5506  in  FIG. 55 ), and the need to maintain previously computed state  5522  and  5524  for particular intervals of time  5510  and  5512 . Thus, the repair-based continuous optimization technique allows for a somewhat continuous recalculation of state functions without full reoptimization in a way that minimizes disturbance and disruption of system control and prevents discontinuities in control. 
   The following is a more formal description of the repair method that represents one embodiment of the present invention. The control optimization problem, solutions for which provide the control state functions that are repaired by the above-described repair method, are given by: 
             min     v   ⁡     (   τ   )         ⁢       ∫     t   1         t   1     +   T       ⁢         ϕ   ~     ⁡     (       x   ⁡     (   t   )       ,     v   ⁡     (   t   )       ,     d   ⁡     (   t   )         )       ⁢           ⁢     ⅆ   t               
subject to constraints:
 
                 x   i     ⁡     (   t   )       -       x   i     ⁡     (     t   1     )       -       ∫     t   1     t     ⁢           f   ~     i     (       x   ⁡     (   τ   )       ,     v   ⁡     (   τ   )         ⁢           )     ⁢     ⅆ   τ           =   0                 g   ⁡     (       x   ⁡     (   t   )       ,     v   ⁡     (   t   )         )       ≥   0               v   ∈   V               t   ∈     [       t   1     ,       t   1     +   T       ]           
where V is a compact subset of R r  and d(t) is a driving function (e.g. deterministic or forecasted demand). {tilde over (φ)}, {{tilde over (f)} i , i=1, . . . , n−1}, and {g j , j=1, . . . , m}are sufficiently smooth. A scalar variable x n (t)≧0 is introduced such that:
 
               x   n     ⁡     (   t   )       =       ∫     t   1     t     ⁢       ϕ   ~     ⁢           ⁢     (       x   ⁢           ⁢     (   τ   )       ,     v   ⁢           ⁢     (   τ   )       ,     d   ⁢           ⁢     (   τ   )         )     ⁢     ⅆ   τ                         x   n     ⁡     (     t   1     )       =   0         
The control optimization problem then becomes:
 
             min     v   ⁡     (   t   )         ⁢       x   n     ⁡     (       t   1     +   T     )             
subject to, for t∈[t 1 , t 1 +T)
 
                 x   i     ⁡     (   t   )       -       x   i     ⁡     (     t   1     )       -       ∫     t   1     t     ⁢         f   i     (       x   ⁡     (   τ   )       ,     v   ⁡     (   τ   )       ,     d   ⁡     (   τ   )         ⁢           )     ⁢     ⅆ   τ           =   0                   g   i     ⁡     (       x   ⁡     (   t   )       ,     v   ⁡     (   t   )         )       ≥   0               v   ∈   V               where   ⁢     :                   f   i     =     {               f   ~     i     ⁡     (       x   ⁡     (   t   )       ,     v   ⁡     (   t   )         )                 ⁢         for   ⁢           ⁢   i     =   1     ,   …   ⁢           ,     n   -   1                           ⁢   ϕ     ~     ⁡     (       x   ⁡     (   t   )       ,     v   ⁡     (   t   )       ,     d   ⁡     (   t   )         )               for   ⁢           ⁢   i     =   n                   
The solution of this problem, v*(t), is called the “optimal control” and the corresponding trajectory, x*(t), is called the optimal trajectory. The above-described control state function (e.g.  4906  and  4916  in  FIG. 49 ) is a combination of x(t) and v(t), t∈[t 1 , T).
 
   A previously computed control state function, with respect to the conventions used in  FIGS. 52–56 , is a combination of:
 
v (t     1     −ΔT)* (t)
 
x (t     1     −ΔT)* (t)
 
The computed differences (e.g.  5602  in  FIG. 56 ), called “the repair,” are defined to be a combination of:
 
δx i   t     1   (t)
 
δv k   t     1   (t)
 
Fencing indicator functions specify whether t is in an interval in which no deviations are allowed in the repaired control state function with respect to each control state variable:
 
               S     x   i       t   1       ⁡     (   t   )       =     {             0             ⁢     if   ⁢           ⁢   t   ⁢           ⁢   is   ⁢           ⁢   in   ⁢           ⁢   the   ⁢           ⁢     “   fenced   ”     ⁢           ⁢   interval   ⁢           ⁢   for   ⁢           ⁢     x   i                 1       otherwise         ⁢     
     ⁢       S     v   k       t   1       ⁡     (   t   )         =     {         0             ⁢     if   ⁢           ⁢   t   ⁢           ⁢   is   ⁢           ⁢   in   ⁢           ⁢   the   ⁢           ⁢     “   fenced   ”     ⁢           ⁢   interval   ⁢           ⁢   for   ⁢           ⁢     v   k                 1       otherwise                     
Assuming a finite number of fenced intervals in the interval [t 1 ,t 1 +T], then for (i=1, . . . , n), where n is the dimension of the model, the incremental state trajectory δx i   t     1   (t) is given by:
 
             δ   ⁢           ⁢       x   i     t   1       ⁡     (   t   )         =         x   i     t   1       ⁡     (     t   1     )       -       x   i       (       t   1     -     Δ   ⁢           ⁢   T       )     *       ⁡     (     t   1     )       -     A   i     t   1       +         S     x   i       t   1       ⁡     (   t   )       ⁢       ∫     t   1     t     ⁢       [             ∑     j   =   1     n     ⁢         ∂     f   i         ∂     x   j         ⁢     (       (       t   1     -     Δ   ⁢           ⁢   t       )     *     )     ⁢       S     x   j       t   1       ⁡     (   τ   )       ⁢   δ   ⁢           ⁢       x   j     t   1       ⁡     (   τ   )                     +       ∑     k   =   1     r     ⁢         ∂     f   i         ∂     v   k         ⁢     (       (       t   1     -     Δ   ⁢           ⁢   t       )     *     )     ⁢       S     v   k       t   1       ⁡     (   τ   )       ⁢   δ   ⁢           ⁢       v   k     t   1       ⁡     (   τ   )                       +       ∑     l   =   1     n     ⁢         ∂     f   i         ∂     d   l         ⁢     (       (       t   1     -     Δ   ⁢           ⁢   t       )     *     )     ⁢   δ   ⁢           ⁢       d   l     t   1       ⁡     (   τ   )                   ]     ⁢           ⁢     ⅆ   τ                   
where the notation (t 1 −Δt)* means that the partial derivatives are evaluated along the reference trajectory.
 
The set {[t s   i,a ,t s   i,b ], s=1, . . . , M i } is the fenced time sub-intervals for state i=1, . . . , n. All of the events A i   e δ(t−t e ) that happen during the time interval [t 1 −ΔT, t 1 ] are approximated with a single event A i   t     1   δ(t−t 1 ) at t 1  as follows:
 
               ∫     t   1   -       t   1   +       ⁢         e   i     ⁡     (   τ   )       ⁢           ⁢     ⅆ   τ         =         ∫     t   1   -       t   1   +       ⁢       A   i     t   1       ⁢     δ   ⁡     (     τ   -     t   1       )       ⁢           ⁢     ⅆ   τ         =     A   i     t   1               
δx i   t     1   (t 1 )=x i   t     1   (t 1 )−x i   (t     1     −ΔT)* (t 1 )−A i   t     1    provides the initial conditions for δx i   t     1   (t)
 
Notice that v (t   −ΔT)* (t) and x (t     1     −ΔT)* (t) are defined on [t 1 −ΔT, t 1 +T−ΔT]. As discussed above, the control state function needs to be extended to the interval [t 1 , t 1 +T]. Assuming that the extension is due to a constant control, and not allowing fencing in the extension:
 
             δ   ⁢           ⁢       x   i     t   1       ⁡     (     t   1     )         =       δ   ⁢           ⁢       x   i     t   1       ⁡     (       t   1     +   T   -     Δ   ⁢           ⁢   T       )         +       ∫       t   1     +   T   -     Δ   ⁢           ⁢   T       t     ⁢       [             ∑     j   =   1     n     ⁢         ∂     f   i       (       t   1     -     Δ   ⁢           ⁢   T       )     *           ∂     x   j         ⁢     (       t   1     +   T   -     Δ   ⁢           ⁢   T       )     ⁢   δ   ⁢           ⁢       x   j     t   1       ⁡     (   τ   )                     +       ∑     k   =   1     n     ⁢         ∂     f   i       (       t   1     -     Δ   ⁢           ⁢   T       )     *           ∂     v   k         ⁢     (       t   1     +   T   -     Δ   ⁢           ⁢   T       )     ⁢   δ   ⁢           ⁢         v   _     k     t   1       ⁡     (   τ   )                       +       ∑     l   =   1     n     ⁢         ∂     f   i       (       t   1     -     Δ   ⁢           ⁢   T       )     *           ∂     d   l         ⁢     (       t   1     +   T   -     Δ   ⁢           ⁢   T       )     ⁢   δ   ⁢           ⁢       d   l     ⁡     (   τ   )                   ]     ⁢           ⁢     ⅆ   τ                   t∈[t   1   +T−ΔT, t   1 +T], i=1, . . . , n 
   where δ  v   t     1   (t)=v(t 1 +T)−v*(t 1 +T−ΔT)* 
   The dynamics on interval [t 1 , t 1 +T] then are given by: 
             δ   ⁢           ⁢       x   i     t   1       ⁡     (   t   )         =     {               x   i     t   1       ⁡     (     t   1     )       -       x   i       (       t   1     -     Δ   ⁢           ⁢   T       )     *       ⁡     (     t   1     )       -     A   i     t   1                       +       S     x   i       t   1       ⁡     (   t   )         ⁢       ∫     t   1     t     ⁢       [             ∑     j   =   1     n     ⁢         ∂     f   i         ∂     x   j         ⁢     (       (       t   1     -     Δ   ⁢           ⁢   t       )     *     )     ⁢       S     x   j       t   1       ⁡     (   τ   )       ⁢   δ   ⁢           ⁢       x   j     t   1       ⁡     (   τ   )                     +       ∑     k   =   1     r     ⁢         ∂     f   i         ∂     v   k         ⁢     (       (       t   1     -     Δ   ⁢           ⁢   t       )     *     )     ⁢       S     v   k       t   1       ⁡     (   τ   )       ⁢   δ   ⁢           ⁢       v   k     t   1       ⁡     (   τ   )                       +       ∑     l   =   1     n     ⁢         ∂     f   i         ∂     d   l         ⁢     (       (       t   1     -     Δ   ⁢           ⁢   t       )     *     )     ⁢   δ   ⁢           ⁢       d   l     t   1       ⁡     (   τ   )                   ]     ⁢           ⁢     ⅆ   τ           ,     t   ∈     [       t   1     ,       t   1     +   T   -     Δ   ⁢           ⁢   T         ]                       ⁢         δ   ⁢           ⁢       x   i     t   1       ⁡     (       t   1     +   T   -     Δ   ⁢           ⁢   T       )         +       ∫       t   1     +   T   -     Δ   ⁢           ⁢   T       t     ⁢       [             ∑     j   =   1     n     ⁢         ∂     f   i       (       t   1     -     Δ   ⁢           ⁢   T       )     *           ∂     x   j         ⁢     (       t   1     +   T   -     Δ   ⁢           ⁢   T       )     ⁢   δ   ⁢           ⁢       x   j     t   1       ⁡     (   τ   )                     +       ∑     k   =   1     n     ⁢         ∂     f   i       (       t   1     -     Δ   ⁢           ⁢   T       )     *           ∂     v   k         ⁢     (       t   1     +   T   -     Δ   ⁢           ⁢   T       )     ⁢   δ   ⁢           ⁢         v   _     k     t   1       ⁡     (   τ   )                       +       ∑     l   =   1     n     ⁢         ∂     f   i       (       t   1     -     Δ   ⁢           ⁢   T       )     *           ∂     d   l         ⁢     (       t   1     +   T   -     Δ   ⁢           ⁢   T       )     ⁢   δ   ⁢           ⁢       d   l     ⁡     (   τ   )                   ]     ⁢           ⁢     ⅆ   τ           ,     t   ∈     [         t   1     +   T     ,       t   1     +   T       ]                         
Expressing min  v(t) x n (t 1 +T) as Φ(x(t 1 +T)), expansion up to the second order gives:
 
             1   2     ⁢       (     δ   ⁢           ⁢     x   ⁡     (       t   1     +   T     )         )     T     ⁢       Φ   xx     ⁡     (     x   ⁡     (       t   1     +   T     )       )       ⁢   δ   ⁢           ⁢     x   ⁡     (       t   1     +   T     )             
The goal of repair is to minimize the above criterion while minimizing the change with respect to a nominal trajectory. The combined criterion is:
 
               1   2     ⁢       (     δ   ⁢           ⁢       x     t   1       ⁡     (       t   1     +   T     )         )     T     ⁢       Φ   xx     ⁡     (     x   ⁡     (       t   1     +   T     )       )       ⁢   δ   ⁢           ⁢       x     t   1       ⁡     (       t   1     +   T     )         +       1   2     ⁢       (     δ   ⁢           ⁢       x     t   1       ⁡     (   t   )         )     T     ⁢   Q   ⁢           ⁢   δ   ⁢           ⁢       x     t   1       ⁡     (   t   )         +       1   2     ⁢       (     δ   ⁢           ⁢       v     t   1       ⁡     (   t   )         )     T     ⁢   R   ⁢           ⁢   δ   ⁢           ⁢       v     t   1       ⁡     (   t   )               
where Q and R are constant positive definite matrices specified by a user that define a balance between satisfaction of the enterprise criterion and minimizing change.
 
The constraints for the repair problem are:
   g   j ( x   (t     1     −ΔT)* ( t )+δ x   t     1   ( t ),  v   (t     1     −ΔT)* ( t )+δ v   t     1   ( t ))≧0, j=1, . . . , m 
When expanded to the first order:
 
               g   j     ⁡     (           x       (       t   1     -     Δ   ⁢           ⁢   T       )     *       ⁡     (   t   )       +     δ   ⁢           ⁢       x     t   1       ⁡     (   t   )           ,         v       (       t   1     -     Δ   ⁢           ⁢   T       )     *       ⁡     (   t   )       +     δ   ⁢           ⁢       v     t   1       ⁡     (   t   )             )       ≈         g   j     ⁡     (         x       (       t   1     -     Δ   ⁢           ⁢   T       )     *       ⁡     (   t   )       ,       v       (       t   1     -     Δ   ⁢           ⁢   T       )     *       ⁡     (   t   )         )       +         ∂     g   i         ∂   x       ⁢     (         x       (       t   1     -     Δ   ⁢           ⁢   T       )     *       ⁡     (   t   )       ,       v       (       t   1     -     Δ   ⁢           ⁢   T       )     *       ⁡     (   t   )         )     ⁢   δ   ⁢           ⁢       x     t   1       ⁡     (   t   )         +         ∂     g   i         ∂   v       ⁢     (         x       (       t   1     -     Δ   ⁢           ⁢   T       )     *       ⁡     (   t   )       ,       v       (       t   1     -     Δ   ⁢           ⁢   T       )     *       ⁡     (   t   )         )     ⁢   δ   ⁢           ⁢       v     t   1       ⁡     (   t   )                           a   .   e   .     ,     j   =   1     ,   …   ⁢           ,   m     ⁢                 
In order to satisfy the constraint inequality, given that
   g   j ( x   (t     1     −ΔT)* ( t ),  v   (t     1     −ΔT)* ( t ))≧0, j=1, . . . , m 
the second and third terms of the expanded constraint must satisfy
 
   
     
       
         
           
             
               
                 
                   
                     
                       ∂ 
                       
                         g 
                         i 
                       
                     
                     
                       ∂ 
                       x 
                     
                   
                   ⁢ 
                   
                     ( 
                     
                       
                         
                           x 
                           
                             
                               ( 
                               
                                 
                                   t 
                                   1 
                                 
                                 - 
                                 
                                   Δ 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   T 
                                 
                               
                               ) 
                             
                             * 
                           
                         
                         ⁡ 
                         
                           ( 
                           t 
                           ) 
                         
                       
                       , 
                       
                         
                           v 
                           
                             
                               ( 
                               
                                 
                                   t 
                                   1 
                                 
                                 - 
                                 
                                   Δ 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   T 
                                 
                               
                               ) 
                             
                             * 
                           
                         
                         ⁡ 
                         
                           ( 
                           t 
                           ) 
                         
                       
                     
                     ) 
                   
                   ⁢ 
                   δ 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     
                       x 
                       
                         t 
                         1 
                       
                     
                     ⁡ 
                     
                       ( 
                       t 
                       ) 
                     
                   
                 
                 + 
                 
                   
                     
                       ∂ 
                       
                         g 
                         i 
                       
                     
                     
                       ∂ 
                       v 
                     
                   
                   ⁢ 
                   
                     ( 
                     
                       
                         
                           x 
                           
                             
                               ( 
                               
                                 
                                   t 
                                   1 
                                 
                                 - 
                                 
                                   Δ 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   T 
                                 
                               
                               ) 
                             
                             * 
                           
                         
                         ⁡ 
                         
                           ( 
                           t 
                           ) 
                         
                       
                       , 
                       
                         
                           v 
                           
                             
                               ( 
                               
                                 
                                   t 
                                   1 
                                 
                                 - 
                                 
                                   Δ 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   T 
                                 
                               
                               ) 
                             
                             * 
                           
                         
                         ⁡ 
                         
                           ( 
                           t 
                           ) 
                         
                       
                     
                     ) 
                   
                   ⁢ 
                   δ 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     
                       v 
                       
                         t 
                         1 
                       
                     
                     ⁡ 
                     
                       ( 
                       t 
                       ) 
                     
                   
                 
               
               ≥ 
               0 
             
             , 
             
               
 
             
             ⁢ 
             
               a 
               . 
               e 
               . 
             
             , 
             
               j 
               = 
               1 
             
             , 
             … 
             ⁢ 
             
                 
             
             , 
             m 
           
           ⁢ 
           
               
           
         
       
     
   
   Meta-Control and Heirarchical Process Construction 
   The discussion, up to this point, has focused on optimization and adjustment of controlling state functions generated by optimization methods at subsequent points in time. However, the present invention encompasses a far more general and powerful method for constructing complex, hierarchical computational systems and models, and for controlling general computation, whether simple or complex. An example of this general method for computational control has been presented above, when additional constraints are added to the original optimization problem in order to produce the meta-level, minimal-time-control-problem problem objective functional  , and the Hamiltonian for the meta-level problem is minimized in order to generate a near-optimal control policy. 
     FIG. 57  illustrates a single level of meta-control of a computational process. In  FIG. 57 , a first computational problem, or computational model,  5702  appears at a lowest level along a meta-level axis  5704 . The computational problem, or model, is discrete, and can be formally described by an iteration {tilde over (x)} k+1 ={tilde over (x)} k +{tilde over (f)}(x k , u k ). As discussed above introducing the continuous-iteration-variable parameterization of the discrete computational problem, or model, the original computational problem can be viewed as continuous, in order to allow for powerful continuous mathematical methods to be applied in order to solve the problem. Thus, as shown in  FIG. 57 , the problem may be viewed as discrete, on one face  5706 , and as continuous, 
                 ⅆ     x   ⁡     (   τ   )           ⅆ     (   τ   )         =     f   ⁡     (       x   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )         ,         
on another face  5708 . The discrete-to-continuous transformation is an important enabling feature of the meta-control method. Using the continuous view of the problem, additional constraints may be added to form a new problem, such as the minimum-time control problem constructed in the above-described optimization problem. An important aspect of this ability to add constraints is that the newly added constraints are not limited or necessarily tied to the original problem. Again, in the above-described embodiment, the initial problem involves finding an optimal control solution for an optimization problem, but the added constraints relate to finding near-optimal control policies for steering the computation towards a time-efficient computation of the original optimization problem. Other types of constraints might be added. For example, one might wish to add meta-level constraints dealing with days of the week when optimizations should be performed, or constraints that steer the computation towards both time and memory-usage efficiencies. The added constraints thus represent a meta-level view superimposed on the original problem. The ability to add new constraints facilitates a hierarchical computational approach to solving problems, similar to the levels of conceptual design used in designing complex software—such as the well-known top-down approach to problem solving.
 
   In  FIG. 57 , addition of constraints is shown as a first part of a transformation represented by the large arrow  5710 . Notice that the entire original problem is carried forward into the new problem both within the constraints for the new problem and within the terms of the new objective function. Thus, the new problem includes the original problem, and solution of the original problem is a subset of the solution of the new problem. Thus, although at the meta-level, one can focus on new constraints and an entirely different view of the original problem, the original problem is encoded in the new problem. In a second part of the transformation, the Hamiltonian is used in order provide an expression for determining a discrete computational method for computing a meta-control policy. The discrete meta-control policy expression  5712  can itself be viewed as a continuous problem  5714  with introduction of a new continuous iteration variable σ. 
     FIG. 58  illustrates multiple levels of meta-control and/or hierarchical construction of a complex computational process. As shown in  FIG. 58 , repeated meta-level transformations may lead from an original, discrete problem  5802  through an arbitrary number of meta-level problems along a meta-level axis  5810 . As described above, each new meta-level allows for addition of new constrains and creation of new problems, while carrying forward the original problem and all lower meta-levels into the higher-level meta levels. As an example, consider a computational modeling problem of modeling a living cell, a current area of much interest in the scientific community. The cell may be viewed hierarchically as a large number of conceptual levels, beginning with a molecular level governed by quantum mechanics and statistical mechanics, and ascending through supermolecular complexes, macromolecular structures and compartments, organelles, genetic control, environmental influences, and myriad other conceptual levels, each viewed in a different way and as being driven by different kinds of forces and principles. However, each level contributes to the whole of the cell, and the effects of molecular-level events and conditions may profoundly influence higher conceptual levels. For example, spontaneous conversion of a particular nucleic-acid base to a different nucleic acid base within a chromosome may lead to an altered gene that, in turn, may lead to an altered gene product that is a component of a crucial organizing structure or temporal control mechanism that, in turn, may severely impact the viability of the entire cell. Moreover, various portions of the problems at each level may be naturally expressed in integer, Boolean, or other types of variables, while other portions of the problems may be defined in terms of continuous variables. 
   The methods of the present invention are applicable to developing a complex computational model for a system such as a living cell or a real-time distributed operating system. The methods discussed above for transforming non-real variables to real-number variables and then to floating point variables provides a way to incorporate the natural expression of different sub-problems into an all-floating-point problem space to which the methods of the present invention may be applied. The ability to define different meta-level views of the cell-modeling problem, while carrying forward the lower-level views, allows for a hierarchical modeling that preserves the full complexity of lower-level views without requiring the lower-level views to be explicitly revisited at higher meta-levels. 
   More formally, assume that an iteration algorithm is provided as follows:
 
 {tilde over (x)}   k+1   ={tilde over (x)}   k   +{tilde over (f)} ( x   k   , u   k )
 
where u k  are control variables and x k  are decision variables.
 
The algorithm converges to
 
               lim     k   -&gt;   ∞       ⁢     x   k       =     x   _           
The specifics of the algorithm are encoded in the function {tilde over (f)} and the parameters u k .
 
In a continualization procedure, the iteration is converted to a controlled differential equation in terms of x(τ) and u(τ) such that
 
                 lim     k   -&gt;   ∞       ⁢     x   ⁡     (   τ   )         ≈     x   _       =       lim     k   -&gt;   ∞       ⁢       x   k     .             
The variable τ is a continuous iteration variable corresponding to the discrete iteration variable k.
 
The differential equation is of the form:
 
               ⅆ     x   ⁡     (   τ   )           ⅆ   τ       =     f   ⁡     (       x   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )             
for all u(·)∈U(·), where U(·) is the set of feasible parameter values.
 
The optimal control problem is formulated as follows:
 
               minimize       u   ⁡     (   r   )       ,   T       ⁢       ∫   0   T     ⁢       ϕ   ⁡     (       x   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )       ⁢     ⅆ   τ           +     Ψ   ⁡     (     x   ⁡     (   τ   )       )                     subject   ⁢           ⁢   to   ⁢           ⁢       ⅆ     x   ⁡     (   τ   )           ⅆ   τ         =     f   ⁡     (       x   ⁡     (   τ   )       ,     u   ⁡     (   τ   )         )                     u   ⁡     (   τ   )       ∈     U   ⁡     (   τ   )             
The functions φ and Ψ are chosen by the algorithm designer in order to achieve meta objectives.
 
Using the minimum principle of Pontryagin, and with the Hamiltonian given by
   H ( x (τ),  u (τ),  p (τ))=φ( x ( r ),  u (τ))+ p (τ) T   {tilde over (f)} ( x (τ),  u (τ)) 
where p(τ) is the costate, the necessary conditions for x*(τ), u*(τ), ∈[0, T) to be optimal are:
 
               ⅆ       x   *     ⁡     (   τ   )           ⅆ   τ       =     (       ∂     H   ⁡     (         x   *     ⁡     (   τ   )       ,       u   *     ⁡     (   τ   )       ,     p   ⁡     (   τ   )         )           ∂   p       )                     ⅆ     p   ⁡     (   τ   )           ⅆ   τ       =     -       (       ∂     H   ⁡     (         x   *     ⁡     (   τ   )       ,       u   *     ⁡     (   τ   )       ,     p   ⁡     (   τ   )         )           ∂   x       )     T                     p   ⁡     (   T   )       =       ∂     Ψ   ⁡     (     x   ⁡     (   T   )       )           ∂   x                       x   *     ⁡     (   0   )       =     x   0                     u   *     ⁡     (   τ   )       ∈     U   ⁡     (   τ   )             
where the Hamiltonian has an absolute minimum as a function of u at u*(r):
   H ( x *(τ),  u *(τ),  p (τ))≦ H ( x *(τ),  u (τ),  p (τ)), ∀ u *(τ)∈ U (τ), τ∈[0 , T ) 
A convergence sequence {û s (τ)} is generated by:
   û   s+1 (τ)= û   s (τ)+ Ŵ ( x (τ), û s (τ),  p   (τ))   
such that
 
                 lim     s   →   ∞       ⁢         u   ⋓     s     ⁡     (   τ   )         =       u   *     ⁡     (   τ   )         ,       when   ⁢           ⁢   x   ⁢           ⁢     (   τ   )       =       x   *     ⁡     (   τ   )         ,     τ   ∈     [     0   ,   T     )             
The above convergence sequence generation expression is analogous to the original iteration algorithm
   {tilde over (x)}   k+1   ={tilde over (x)}   k   +{tilde over (f)} ( x   k   , u   k ) 
where û s (τ) is analogous to {tilde over (x)} k .
 
Using the continualization procedure described above with σ being a continuous parameterization of s, one obtains:
 
                       ⁢     ∂     v   ⁡     (     τ   ,   σ     )             ∂   σ       =     W   ⁡     (       x   ⁡     (   τ   )       ,     v   ⁡     (   τ   )       ,     p   ⁡     (   τ   )         )             
where v(τ, σ) is a continulaized version of û s (τ) and W(x(τ), v(τ), p(τ)) is a continualized version of Ŵ(x(τ), û s (τ), p(τ)).
 
   To summarize: 
   
     
       
         
           
             
               ∂ 
               
                 v 
                 ⁡ 
                 
                   ( 
                   
                     ɛ 
                     , 
                     G 
                   
                   ) 
                 
               
             
             
               ∂ 
               G 
             
           
           = 
           
             w 
             ⁡ 
             
               ( 
               
                 
                   x 
                   ⁡ 
                   
                     ( 
                     τ 
                     ) 
                   
                 
                 , 
                 
                   v 
                   ⁡ 
                   
                     ( 
                     τ 
                     ) 
                   
                 
                 , 
                 
                   p 
                   ⁡ 
                   
                     ( 
                     τ 
                     ) 
                   
                 
               
               ) 
             
           
         
       
     
     
       
         
           
             
               ∂ 
               
                 p 
                 ⁡ 
                 
                   ( 
                   τ 
                   ) 
                 
               
             
             
               ∂ 
               τ 
             
           
           = 
           
             - 
             
               
                 ( 
                 
                   
                     ∂ 
                     
                       H 
                       ⁡ 
                       
                         ( 
                         
                           
                             x 
                             ⁡ 
                             
                               ( 
                               τ 
                               ) 
                             
                           
                           , 
                           
                             u 
                             ⁡ 
                             
                               ( 
                               τ 
                               ) 
                             
                           
                           , 
                           
                             p 
                             ⁡ 
                             
                               ( 
                               τ 
                               ) 
                             
                           
                         
                         ) 
                       
                     
                   
                   
                     ∂ 
                     x 
                   
                 
                 ) 
               
               T 
             
           
         
       
     
     
       
         
           
             
               ∂ 
               
                 v 
                 ⁡ 
                 
                   ( 
                   
                     τ 
                     , 
                     σ 
                   
                   ) 
                 
               
             
             
               ∂ 
               σ 
             
           
           = 
           
             W 
             ⁡ 
             
               ( 
               
                 
                   x 
                   ⁡ 
                   
                     ( 
                     τ 
                     ) 
                   
                 
                 , 
                 
                   v 
                   ⁡ 
                   
                     ( 
                     τ 
                     ) 
                   
                 
                 , 
                 
                   p 
                   ⁡ 
                   
                     ( 
                     τ 
                     ) 
                   
                 
               
               ) 
             
           
         
       
     
     
       
         
           
             with 
             ⁢ 
             
               
                   
               
               ⁢ 
               
                   
               
             
             ⁢ 
             
               p 
               ⁡ 
               
                 ( 
                 T 
                 ) 
               
             
           
           = 
           
             
               
                 ∂ 
                 
                   Ψ 
                   ⁡ 
                   
                     ( 
                     
                       x 
                       ⁡ 
                       
                         ( 
                         T 
                         ) 
                       
                     
                     ) 
                   
                 
               
               
                 ∂ 
                 x 
               
             
             ⁢ 
             
                 
             
             ⁢ 
             and 
             ⁢ 
             
                 
             
             ⁢ 
             
               x 
               ⁡ 
               
                 ( 
                 0 
                 ) 
               
             
             ⁢ 
             
                 
             
             ⁢ 
             
               given 
               . 
             
           
         
       
     
   
   Although the present invention has been described in terms of a particular embodiment, it is not intended that the invention be limited to this embodiment. Modifications within the spirit of the invention will be apparent to those skilled in the art. For example, an almost limitless number of different implementations using different operating-system platforms, programming languages, data structures, control structures, and other such variations are possible. 
   The foregoing description, for purposes of explanation, used specific nomenclature to provide a thorough understanding of the invention. However, it will be apparent to one skilled in the art that the specific details are not required in order to practice the invention. In other instances, well-known circuits and devices are shown in block diagram form in order to avoid unnecessary distraction from the underlying invention. Thus, the foregoing descriptions of specific embodiments of the present invention are presented for purposes of illustration and description; they are not intended to be exhaustive or to limit the invention to the precise forms disclosed, obviously many modifications and variations are possible in view of the above teachings. The embodiments were chosen and described in order to best explain the principles of the invention and its practical applications and to thereby enable others skilled in the art to best utilize the invention and various embodiments with various modifications as are suited to the particular use contemplated. It is intended that the scope of the invention be defined by the following claims and their equivalents: