Stable adaptive control using critic designs

[0002] The present application is related to and claims priority to U.S. provisional application serial No. 60/101,680, entitled “Stable Adaptive Control Using New Critic Designs,” filed on Sept. 24, 1998, naming Paul J. Werbos as the sole inventor. (The contents of that application are incorporated herein by reference.) The present application is related to U.S. Pat. No. 6,424,056 and co-pending U.S. applications Ser. No. 09/147,338, filed May 10, 1999, and Ser. No. 60/127,952, filed Apr. 6, 1999, all naming Paul J. Werbos as sole inventor. The present application is related to a U.S. Provisional application serial No. 60/152,167 filed Sep. 4. 1999, entitled “NEURAL NETWORKS FOR INTELLIGENT CONTROL,” also naming Paul J. Werbos as sole inventor. The contents of those co-pending applications are incorporated herein by reference.

[0003] 1. Field of the Invention

[0004] The present invention is directed to a method of providing stable control of a system (e.g., an electrical power grid, a factory, and a financial prediction system) using a neural network-based design, and more particularly, to a neural network-based control system and method using critics.

[0005] 2. Discussion of the Background

[0006] The field of intelligent control is an expansive area that has attempted to solve many complex problems. Significant research has previously been performed in the areas of classical control theory [1-6] and biologically-inspired intelligent control [7-15]. As an outgrowth of that and other research, researchers have attempted to solve the Generalized Moving Target (GMT) problem [16,17], which is defined as follows: for a function E(v, w), find two vectors v and w such that the following two conditions are met: (1) E is minimized with respect to v for w fixed; (2)w=v.

[0007] Other research includes two additional classes of control design. The first class is “adaptive control,” in the tradition of Narendra, which includes both linear (multiple-input multiple-output, MIMO) designs [1] and nonlinear or neural extensions [18,19]. The second class is learning-based Approximate Dynamic Programming (ADP) [7-15], which has sometimes been presented as a form of “reinforcement learning” [20-22], sometimes called “adaptive critics” [23] and sometimes called “neuro-dynamic programming” [24].

[0008] Previous forms of adaptive control discussed by Narendra, had difficulty in ensuring stability even in the linear-quadratic case. Those designs appear very similar to a particular ADP design, HDP+BAC[8,22,25], that was developed back in 1971, and reported in Werbos' Ph.D. proposal[26]. In fact, when HDP+BAC is applied to the linear-quadratic case, it reduces down to a standard indirect adaptive control (IAC) design, except that there is a provision for adapting one extra matrix which, in principle, should permit much stronger stability guarantees. Roughly speaking, ordinary LAC designs are based on the minimization of tracking error, which is sometimes just: $\begin{array}{cc}\sum _i{\left(x_i^{\left(\mathrm{ref}\right)}\left(t+1\right)-x_i\left(t+1\right)\right)}^2={\underset{\_}{e}\left(t+1\right)}^T\underset{\_}{e}\left(t+1\right),& \left(1\right)\end{array}$

[0009] where the vector x^(ref)(t) represents the desired state of the plant or external environment at time t, where x(t) represents the actual state, and e(t) represents the gap or error between the two. In the same situation, HDP+BAC reduces to the same overall design, except that:

[0010] Ĵ(t+1)=e(t+1)^TCe(t+1),  (2)

[0011] is minimized, where C is a matrix of “Critic weights” to be adapted by HDP. In the nonlinear case, equation (2) is replaced by an artificial neural network (ANN), in order to allow Ĵ to approximate any nonlinear function. (In HDP and GDHP, it is possible to supplement the training of the Action network, to include some stochastic search (particularly in the offline mode) to help keep the system out of local minima [46]. Barto, Thrun and others have discussed other approaches to adding noise to the Action network or controller [8,21].)

[0012] In theory, if the critic weights C or the Critic network converged to the right values (the values which satisfy the Hamilton-Jacobi-Bellman equation [4,5,24,28]), then the function Ĵ would serve as a Liapunov function guaranteed to stabilize the overall system, if the system is controllable. However, none of the established methods for adapting a Critic possess quadratic unconditional stability, even in the linear deterministic case. The variations of these methods based on the Galerkin approach to solving differential equations do possess unconditional stability in that case, but they converge to the wrong weights almost always in the linear stochastic case.

[0013] The notation herein generally tracks the sources upon which it draws. However two biased selections have been made herein. Herein, “C”, “b” and W are used for the estimated values of certain sets of parameters or weights, and “C*” or “b*” are used for their true values. On the contrary, for the functions J and λ, the usual three-fold convention of using the “J” and “λ” for the true values, “J*” and “λ*” for target values used in adaptation, and “Ĵ” and “{umlaut over ({circumflex over (e)})}” for the estimated values are used. The state vector is identified as “∂_tx” in linear adaptive control, rather than the usual x-dot. The notation “∂_t” has long been part of the standard notation in physics.

[0014] In principle, the ideal feedback control system should try to optimize some mix of stability and performance, in the face of three general types of uncertainty:

[0015] (1) High-bandwidth random disturbances, which are usually represented as a stochastic process, based on random noise [4,5,26,29], but are sometimes represented as bounded disturbances of unknown structure [1];

[0016] (2) Drifting values (and occasional abrupt changes) in familiar process parameters such as friction, viscosity and mass, often due to the aging of a plant or the changes in general environmental conditions:

[0017] (3) Uncertainty about the fundamental structure of the plant or environment (sometimes due to catastrophic events, like a wing being shot off of an airplane), and shifts of parameters in ways that could not be anticipated or simulated even as possibilities at the time when the controller is developed.

[0018] This three-fold distinction is difficult to formalize, in mathematical terms, but it has great practical importance. Roughly speaking, the ability to respond to the first type of disturbance or uncertainty may be called “stochastic feedback control.” The ability to respond to the second type may be called “adaptation,” and the ability to respond to the third type is called “learning.” The practical tradeoffs here are discussed at some length in the introductory review in [7], and in other papers cited in [7].

[0019] “Adaptive control” [1-3] has often been viewed as a tool for addressing the second type of uncertainty—uncertainty about drifting plant parameters. The most classical designs in adaptive control are intended to control plants x(t) governed by the equations:

[0020] ∂_t=A  (3)

[0021] or:

[0022] (t+1)=A(t)+B(t),  (4)

[0023] where u is a vector of controls, where A and B are unknown matrices representing the parameters of the plant, and where ∂_trepresents differentiation with respect to time. In the

[0024] simplest case, the state vector x is directly observable. The key idea is to develop control designs which estimate the matrices A and B, explicitly or implicitly, as the process unrolls in real time, and converge “on the fly” to a good control strategy:

[0025] (t)=K(t)(t),  (5)

[0026] despite the ignorance about A and B. In the general case, it is assumed that x(t) is not directly observable. Instead, a vector v governed by:

[0027] (t)=H(t)  (6)

[0028] is observed. Roughly speaking, this requires the use of a more complex control strategy like:

[0029] (t)=K_1,1(t−1)+K_1,2(t−2)+ . . . +K_1,k(t−k)+K_2,1(t−1)+ . . . +K_2,k(t−k),  (7)

[0030] for some integer k. (See [1, p.411].) There exists a huge body of stability theorems for adaptive control, both in the original linear versions [1-3] and in a variety of nonlinear and neural extensions (e.g., [18,19]).

[0031] In practice, however, all of these stability theorems require very strong assumptions about the plant or environment to be controlled. Ordinary adaptive control and neural adaptive control have often exhibited problems with stability and with slow response to transient disturbances, particularly in real-world plants containing delays and deadzones and reversal phenomena, etc.

[0032] Because of these problems, the most powerful approaches available today to cope with type-two uncertainty in real engineering applications are:

[0033] (1) Adaptive-predictive control [2], in the form which explicitly minimizes tracking error over multiple time periods ahead into the future. This involves a substantial increase in computational complexity, for computations that must be performed in real time.

[0034] (2) A two step design process, performed offline before the controller is actually used on the real commercial plant. In the first step, one designs a controller containing very strong feedback loops and/or observers, able to estimate (explicitly or implicitly) the specific plant parameters that are expected to drift. These loops take over the role which the estimates of “A” and “B” would play in an adaptive control approach. In the second step, one exploits prior knowledge about the specific plant, including knowledge about its uncertainties, in order to tune or train this controller. The controller is tuned or trained for optimal performance over multiple time periods, in some kind of offline training process.

[0035] Optimization may be done by using (1) an n-period lookahead method, analogous to model-predictive control [2.8], or (2) dynamic programming. The second approach often has very good response to transient disturbances, because of how it exploits prior knowledge about the plant. The first approach, however, appears to have had very limited success. Astrom's work in general [2] is well-respected for the many industrial applications it has led to.

[0036] Athans and Baras [6] have been very strong spokesmen for the second approach, within the nonlinear robust control community. Athans, in particular, has worked with many aerospace applications, where freezing the dynamics of a controller prior to actual flight can simplify the process of obtaining government approval. Robust control theorists have shown that type-one disturbances convert the problem of controller design into a problem of stochastic optimization. (This is true no matter how the disturbances are represented—either as stochastic noise or as “worst case noise” in the spirit of Tsypkin.) Nonlinear stochastic optimization problems require the use of dynamic programming, rather than a simple n-step lookahead. In order to find numerical solutions to such stochastic optimization problems, one can simply use learning-based ADP, applied in an offline learning mode.

[0037] In the neural network field, the second approach was first proposed in 1990 [30] and called “learning offline to be adaptive online.” In one embodiment, a Time-Lagged Recurrent Network (TLRN) is used as a controller [8,ch10; 26.ch.8] with feedback loops in the control network. Additional inputs can be fed to the controller, by including the intermediate outputs of TLRNs trained to perform system identification of the plant. Such intermediate outputs serve as neural observers. Learning offline to be adaptive online is the foundation of Ford's “multistreaming approach” [7,31], which underlies the most successful and impressive success of neural networks to date in real-world control applications. Most of the Ford work has used n-step lookahead optimization, using backpropagation through time (BTT) to minimize the computational costs. With BTT and special chips, Fold expects to be able to extend this work to include adaptation on the fly, so as to upgrade the initial controller developed offline. BTT [26] can be used to reduce the cost of computing exact derivatives through any feedforward differentiable system, not just neural networks. (For systems which are not feedforward, see the review in [32]. See www.nd.com for some tools to implement BTT and various special cases of TLRN.)

[0038] Many researchers have been overwhelmed by the sheer complexity of biological brains. Particularly in the cerebral cortex of mammals, there is a bewildering array of cells that perform a wide variety of tasks: Many people have despaired of finding any universal principles of information processing at work there, which could be understood in mathematical terms. However, studies of mass action in the cerebral cortex have demonstrated that cells in one part of the cortex can learn to take over the functions of other cells, when there is a need to do so (and when the required data inputs are not cut off). (See the work of Lashley, Pribram and Freeman and [9-11,26].) However, the problem lies in developing the right kind of mathematics. Since the overall function of a biological brain is to compute actions or decisions, learning-based intelligent control may someday provide the required mathematics. (See [9-11] for more specific links between ADP designs and the brain.)

[0039] Nevertheless, the role of real-time learning should not be overstated, even in biology. Since the time of Freud at least, it has been well known that organisms remember past experiences, and adapt their general models or expectations about their environment based on some combination of present experience and memories. Back in 1977 [25], a simple learning scheme was proposed based on the more general idea of “syncretism,” an interpretation of Freud's vision of the interplay between memory and generalization. Syncretism promises substantial practical benefits for neural network adaptation [34: 8,ch.3]. Early successes of “memory-based learning” [35]'support this general approach. More recently, McClelland and others have promulgated theories of the interplay between memory and generalization that also embody a form of this idea [36], although the inventor has previously predicted that this interplay mainly occurs within the cerebral cortex and the hippocampus, rather than between them, as suggested by McClelland. Recent arguments by Pribram,. Alkon and others tend to support the former view [36]. But in any case, the psychological evidence for the phenomenon as such appears very convincing. Within the world of artificial intelligence, researchers such as Rosalind Picard of MIT have developed learning methods that also embody this general approach.

[0040] In an ideal world, in the linear deterministic case, a universal adaptive controller would exist with the following properties. First, the design and the adaptation rule would produce the output the vector u(t) at each time t, based only on knowledge about u(τ) and v(τ) at previous times τ<t. It would not require knowledge about A, B or H, except for the dimensionality of u and v (and perhaps x). It should be guaranteed to send tracking error to zero, asymptotically, at an exponential rate, for all combinations of A, B and H and reference model for which this would be possible with a fixed controller, (The “reference model” is the process—generally considered linear here—which outputs the desired observed state v*(t).) In other words, the ideal adaptive controller should possess the same kinds of strong stability guarantees that exist for more familiar kinds of control, such as fixed-stricture LQG optimal control [4,5].

[0041] Simplistic adaptive control designs available today fill far short of this ideal. (Certain hybrid designs of Astrom[2] are more complex and are discussed later.) The bulk of the existing theory for linear adaptive control focuses on the Single-Input Single-Output (SISO) special case—the case where y and u are scalars. Yet even in that case, stability is normally guaranteed only under three restrictions (in addition to restrictions on the reference model). Those restrictions are: (1) information must be available to choose the integer “k” in equation 7; (2) the plant must be “minimum phase”; and (3) the sign of the high-frequency gain must be known Narendra and Annaswamy go on to state [1,p.359]: “ . . . it was realized that these assumptions were too restrictive, since in practice they are violated by most plants, even under fairly benign conditions.”

[0042] The problem with the sign of the gain may be illustrated by a simplified, commonsense example. A fisherman has complete control over the big lake in which he fishes. The “reference model” is simply a certain quota (in weight) of fish that must be caught every year. One month the quota is not met The lake is usually capable of producing enough fish to meet the quota, in the long term, but an unknown disturbance has caused the fish population to go through a minor dieback. Following the usual policy of minimizing tracking error at time t+1, the fisherman increases his catch to arrive back at the quota. However, by fishing at a higher level. the fish population is reduced still further. This dips into the smaller fish that could have provided more growth in the future. Thus to meet the quota in the following month, even more fish are caught. In the end, the entire fish population dies off. That is, things go nonlinear, resulting in a catastrophic instability. This example leads to a host of further implications. But the key point is that policies that reduce tracking error (or maximize utility) in the short term may have undesirable or even catastrophic effects in the more distant future. Knowing the signs of the gains is crucial to avoiding those kinds of problems in adaptive control designs that do not look beyond time t+1. (The “minimum phase” assumption has similar, related pervasive implications, discussed at length by many researchers, including Widrow [37].)

[0043] Narendra and Annaswamy discuss methods that avoid the need to know the sign of the high-frequency gain, for the SISO case. However, those methods use a very special-purpose trick, the Nussbaum gain, which does not appear to carry over to the general, MIMO case in a serious way. The requirements for prior knowledge in the MIMO case [1, ch. 10] are far more complex and demanding than in the SISO case. After all, in the SISO case, there are only two possibilities for the sign of the high-frequency gain, which is a scalar—plus or minus. In the MIMO case, there are an infinite number of possible directions or modes of instability, and the requirements are very complex. Narendra has recently found another way to avoid having to know the sign of the gain, by using multiple models [38,39] in the SISO case. However, the minimum phase assumption is still required, and a universal MIMO controller developed in this way would be extremely complicated, if possible.

[0044] Difficulties with these restrictions probably explain why Lyle Ungar, of the University of Pennsylvania, found unstable results when he tested the most advanced direct, indirect and hybrid direct-indirect adaptive control designs on the bioreactor benchmark test problem given in [22]. There may even be an analogy between the harvesting of cells from a bioreactor and the harvesting of fish from a lake. This same problem was later solved directly and efficiently, both by neural model-predictive-control (based on BTT) and by ADP methods, in papers from Texas Tech [14] and from Ford Motor Company [34].

[0045] Narendra and Annaswamy [1] discuss both direct and indirect adaptive control designs for the linear case. However, in, 1990, Narendra [22,p. 135] stated that until “direct control methods are developed, adaptive control of nonlinear dynamical systems has to be carried out using indirect control methods.” In 1992, he stated [8,p.168] that “unless further assumptions concerning the input-output characteristics of the plant are made, direct adaptive control is not possible.”

[0046] FIG. 1 is essentially identical to the nonlinear adaptive control design discussed by Narendra in [22,p. 135], [8,p.168], [41,p. 166] and elsewhere. (In the neural network context, this application uses the capital letter, X, rather than v, to indicate the vector of observables. R is used to indicate the estimated state vector (or “representation of reality,” usually based on recurrent neurons). Moreover, the phrase “Model network” is used to describe what Narendra calls the identification network, Nⁱ. The phrase term “Action network” is used to describe what he calls the controller network, N^c. Both Narendra and this application assume that the Model network may be adapted in real time, by an adaptation rule independent from what is used to adapt the Action network. However, in his flowcharts, Narendra adds a few arrows (labeled “e_i”) to give a hint of how the Model network might be adapted.

[0047] Numerous implementations of the design in FIG. 1 have been reported in journal articles all over the world. Two points must, however, be known:

[0048] (1) How to adapt the Model network. Numerous ways [8,32] of doing this have been suggested with varying degrees of robustness. For linear adaptive control, simple least mean squares (LMS) learning may be adequate. (For example, see [1,p.402], [42].)

[0049] (2) How to adapt the Action network. Narendra proposes the weights W_ijin the Action network be adapted in proportion to the derivatives of tracking error with respect to the weights. This may be written: $\begin{array}{cc}{W}_{\mathrm{ij}}\left(t+1\right)=W_{\mathrm{ij}}\left(t\right)-\alpha \frac{\partial }{\partial W_{\mathrm{ij}}}\left(E\left(t+1\right)\right),& \left(8\right)\end{array}$

[0050] where:

[0051] E(t+1)=((t+1))²=(*(t+1)−(t+1))²  (9)

[0052] is the tracking error at time t+1 and α is some arbitrary (small, positive) learning rate. Narendra describes calculating the derivatives in equation 8 by “backpropagating through the Model network.” How to do this is described at length both in Narendra's work and—for a broader class of possible network structures (neural or nonneural)—in other references. [8.26]. The broken arrows in FIGS. 1 and 2 represent the backwards calculations used to obtain the required derivatives at minimum computational cost. Additional literature discusses (1) alternative ways to choose the learning rate and (2) alternative gradient-based learning rules.

[0053] FIG. 1 has many limitations, as discussed herein and in other references [41]. Applied to the case of a fully observable linear plant, equation 8 reduces to:

[0054] K(t+1)=K(t)−αB^T(t+1)^T(t),  (10)

[0055] where the weights in the Action network are now just the matrix K of equation 5. “Backpropagating through a network” (like the Model network, which reduces to equation 4 in the linear deterministic case) is simply a low-cost way to multiply a vector by the transpose of the Jacobian of that network; it reduces costs by exploiting the internal structure of the network, and working backwards.

[0056] FIG. 2 is very similar to FIG. 1, in a mechanical sense, but it has vastly different properties. In FIG. 2, the terms “HDP+BAC” refer to Heuristic Dynamic Programming (HDP) and the Backpropagated Adaptive Critic (BAC). HDP is a technique for adapting the Critic network, developed in 1968-1972. BAC is the method used to adapt the Action network. HDP+BAC is an architecture for optimal control, which is structurally almost identical to IAC. HDP+BAC can be used in a variety of configurations, involving a combination of real-time learning, offline learning, prior information, etc. In order to implement a full, real-time learning version of HDP+BAC, three networks must be adapted (i.e., the Critic, Model and Action networks must be adapted).

[0057] Note that the design here is identical to the IAC design in FIG. 1, except that tracking the error (E) has been replaced by the Critic network. The Model network can be adapted by using any of the various methods previously proposed and used for that task. The Action network can be adapted by replacing equation 8 by: $\begin{array}{cc}{W}_{\mathrm{ij}}\left(t+1\right)=W_{\mathrm{ij}}\left(t\right)-\alpha \frac{\partial }{\partial W_{\mathrm{ij}}}\left(\hat{J}\left(t+1\right)\right),& \left(11\right)\end{array}$

[0058] where the derivatives are again calculated by backpropagation, in exactly the same way. (To initialize the backpropagation, however, the derivatives of Ĵ are first calculated with respect to its inputs, R(t+1,). Those derivatives can be calculated quite easily by backpropagating through the Critic network. In the deterministic case, the distinction between R and {circumflex over (R)} is not so important as it is in the stochastic case [8].) Strictly speaking, FIG. 2 represents a special case of HDP+BAC applicable to the tracking problem. In other words, to implement FIG. 2, the same computer code can be used to implement FIG. 1, except that “E” is replaced by “Ĵ.”

[0059] HDP attempts to adapt the Critic network in such a way that its output, Ĵ, converges to a good approximate solution to the Bellman equation of dynamic programming. This equation may be written fairly generally [8] as: $\begin{array}{cc}J\left(\underset{\_}{R}\left(t\right)\right)=\underset{\underset{\_}{u}\left(t\right)}{\mathrm{Max}}\left\{U\left(\underset{\_}{R}\left(t\right),\underset{\_}{u}\left(t\right)\right)+\left(1/\left(1+r\right)\right)⟨J\left(\underset{\_}{R}\left(t+1\right)\right)⟩\right\}-U_0,& \left(12\right)\end{array}$

[0060] where U is the utility or cost function to be maximized or minimized in the long term, where r is an interest rate parameter used to discount the value of fixture utility, where the angle brackets < > indicate expectation value, and where U₀is a parameter introduced by Howard [28] to extend dynamic programming to the case of an infinite time horizon with r=0. (Of course, for a minimization task, “Max” is replaced by “Min” in equation 12.) When this method is applied to pure tracking problems, as in classical adaptive control, U is chosen to be the tracking error E, and the reference model is treated as a fixed augmentation of the Model network.

[0061] Karl Astrom has perhaps been the world's number one leader in stimulating real-world applications of adaptive control. He has done this in part by developing sophisticated hybrid designs, in order to meet the demanding requirements of these applications. For example, he played a pioneering role in building a link between linear/quadratic optimization and adaptive control. In [8,ch.2], he and McAvoy give a broad overview of efforts by themselves and others to address the general subject of intelligent control. The latest edition of his classic text on adaptive control [45] contains many references to new work in these directions.

[0062] The original, classic text by Astrom and Wittenmark [2] discusses three types of adaptive control which overcome many of the restrictive assumptions discussed by Narendra: (1) multiperiod adaptive-predictive control; (2) adaptive pole placement; and (3) linear-quadratic Self-Tuning Regulators (STRs). (Astrom also has previously provided insight to the minimum phase problem. He pointed towards the adaptive pole placement and LQG/STR approach as the most important ways to overcome this problem.)

[0063] The future plans by Ford Motor Company to consider using multiperiod optimization in real time, along with real-time system identification, could be viewed as a possible future application of multi-period adaptive-predictive control. There are many other links, to the most successful applications of artificial neural networks [7]. Methods of this family are of great practical utility, but they (and adaptive pole placement methods) are not “true” adaptive control designs as discussed herein. They are also not plausible as a theory of how the brain actually implements intelligent control. Also, it is hard to imagine a nonlinear stochastic generalization of that design. The discussion in [2] describes the method as an SISO method.

[0064] On the other hand, the discussion of linear-quadratic STRs [2] has many close links to the present invention. One version, based on spectral factorization, does not appear relevant to the goals described herein. But the other version—Indirect STR based on Riccati equation—is extremely close. This version is called Algorithm 5.7 in [2], and Algorithm 4.4 in [45]. The new edition of the text is almost identical to the old version in this section. The matrix S(t) in Astrom's equation 5.47 resemble the matrix of Critic weights C discussed herein. The methods used by Astrom and others to update the estimate of S(t) could be viewed as specialized Critic adaptation methods, for use in the linear/quadratic case. Unfortunately, the most reliable update methods—based on solving a Riccati equation or performing a spectral factorization—are not true adaptive control methods as defined above. Astrom hints that there have been serious stability problems when other update methods have been attempted. However, he also hints that D. W. Clarke of Oxford and M. Karny of Prague have been more central to this work in recent years. Karny has certainly made significant contributions to the larger area of intelligent systems [46].

[0065] Landelius' method of implementing HDP was motivated in part by the most recent work in this literature [47], which he also cites and discusses.

[0066] As an another approach, Liapunov stability theory in general has influenced huge sections of control theory, physics, and many other disciplines. More narrowly, within the disciplines of control theory and robotics, many researchers have tried to stabilize complex systems by first deriving Liapunov functions for those systems. In some cases, the Liapunov functions have been derived analytically by solving the multi-period optimization problem in an analytic fashion.

[0067] Having derived an application-specific Liapunov function, one can then use the design of FIGS. 1 and 2, in which the Liapunov function replaces the square tracking error or J. Theoretically, by replacing the square tracking error with some other pre-specified error measure, new stability properties and new restrictions on the plant are created.

[0068] The Liapunov function approach has been particularly useful in the field of robotics [48], where many robot arms obey nonlinear but rigid dynamics. Robert Sanner of the University of Maryland and John Doyle of CalTech are often mentioned in discussions about the successes of this approach. Nevertheless, none of these analytically derived fixed Liapunov functions can be regarded as a universal adaptive controller, as discussed above, any more than IAC itself can.

[0069] From a practical point of view, it becomes increasingly difficult to derive such Liapunov functions analytically, as the complexity of the nonlinear systems (e.g., elastic, light-weight and flexible robot arms) increases. (See [7] for an explanation of some of Hirzinger's success in working with such robots. Some of the unpublished work of Fukuda has been even more useful.) The difficulties here are analogous to the difficulty of trying to solve simple algebraic equations analytically. As the order of equations increases, eventually a point is reached where closed-form analytic methods simply cannot provide the solution.

[0070] Hundreds of papers have been published by now on various forms of approximate dynamic programming (ADP), adaptive critics and/or reinforcement learning. However, only three of the established, standard methods for adapting Critic networks are directly relevant to the goals herein: (1) Heuristic Dynamic Programming (HDP); (2) Dual Heuristic Programming (DHP); and (3) Globalized DHP (GDHP). The classic “Temporal Difference” (TD) method is essentially a special case of HDP. The recent “two-sample” method [24,50,51] is also relevant, but is discussed later.

[0071] As discussed above, the present invention focuses on achieving stability for methods of adapting Critic networks, not for entire adaptive critic control systems. Accordingly, a complete review of the larger question of how to adapt a complete ADP control system is not provided herein. For the concurrent adaptation of Model networks, Action networks and Critic networks, and for practical experience, see [8,14,15,33,52-56].

[0072] The only existing, working control designs which meet criteria/definition for “Model-Based Adaptive Critics” (MBAC) or “brain-like intelligent control” [33] described herein are those based on FIG. 2, or the equivalent of FIG. 2 for DHP, GDHP and their variations. MBAC has been successfully implemented by at least three companies (Ford Motor Co., Accurate Automation Corp., Scientific Cybernetics Inc.), four professors working with graduate students (Wunsch, S. Balakrishnan, Lendaris, W.Tang), and three other individuals (Jameson, Landelius and Otwell). Since the first early implementations in 1993, MBAC has outperformed other modern control and neurocontrol methods in a variety of difficult simulated problems, ranging from missile interception (Balakrishnan) to preventing cars from skidding when driving over unexpected patches of ice (Lendaris). The one physical implementation published to date was also highly successful (Wunsch's solution of Zadeh's “fuzzy ball and beam” challenge). Balakrishnan has performed work on a larger physical implementation (a cantilever plate). The one alternative method which still seems competitive with MBAC, in terms of performance, in these tests, is neural model-predictive control based on BTT.

[0073] Roughly speaking, the adaptive critic field of research is a single new research field (albeit still divided into factions) which emerged around 1988-1990 through the unification of several previously separate strands of research. The origins of the field up to that time are summarized in FIG. 3. Within the field itself, the terns “reinforcement learning,” “adaptive critics.” “approximate dynamic programming” and “neurodynamic programming” are normally viewed as approximate synonyms. Nevertheless, the choice of terms also reflects different goals for research within the field and different patterns of interest in related research outside of the field. FIG. 3 only represents the flow of key ideas in this topic prior to 1988-1990: it does not represent patterns of personal association, ideas in other areas, or important recent ideas that will be mentioned later.

[0074] The psychologist B. F. Skinner is well-known for the idea that rewards and punishments (“primary reinforcement signals”) determine the behavior of animals including humans. Many of the ADP designs used today in engineering or computer science are also being used as empirical models of animal behavior. Harry Klopf was a major pioneer of that tradition. Klopf; in turn, was responsible for recruiting and funding Barto to explore this strand of research. Barto has continued this tradition through ongoing collaborations with psychologists and neuroscientists.

[0075] Widrow never pursued these kinds of connections, but his seminal 1973 paper on adaptive critics [23] was clearly influenced by Skinner's notion of reward and punishment. Skinner's notion of a “secondary reinforcement system” can be viewed as one way of talking about a Critic network.

[0076] Von Neumann and Morgenstern [57] invented the concept of cardinal utility function, which underlies the ADP approach. This work also made possible many new directions in economic theory, ranging from game theory to decision analysis and Bayesian utilitarianism. Working with Richard Bellman, Von Neumann also helped to inspire the development of dynamic programming. In criticizing the neuron models of McCulloch and Pitts, Von Neumann [58,p.451] made two major points: (1) that neuronal signals may be more accurately modeled as continuous variables rather than binary signals in many cases; (2) that study of mechanisms like learning and memory are more important than efforts to understand how particular static, learned functions can be represented in fixed neural networks. The successful revival of the neural network field in the past decade was based, in large part, on researchers finally embracing these two crucial insights.

[0077] On a philosophical level, Skinner and Von Neumann were extremely far apart. Skinner's words about molding and controlling human beings certainly gave encouragement (even if unintentionally) to a variety of ideologies in the spirit of fascism and Communism. The revival of the neural network field in the 1980's was motivated in part by a rejection of Skinner's approach and the resulting search for new paradigms in psychology [59]. Von Neumann, on the other hand, encouraged a higher level of respect for human intelligence, autonomous human decision-making and human potential, from the very beginning. Despite these differences, tie most basic views and insights of both men have been fully incorporated into the ongoing research in this field. Different researchers have different attitudes about the deeper philosophical implications [12], even as they use the same mathematics.

[0078] The concept of reinforcement learning as a pathway to the construction of intelligent systems has often been credited to the great pioneers of artificial intelligence (AI)—Newell, Shaw and Simon—and to Marvin Minsky [60]. They proposed the development of machines that learn over time to maximize some measure of reward or reinforcement. They proposed that such a machine, simply by learning through experience, could gradually develop higher-order intelligence as a kind of emergent phenomenon. The earliest attempts to implement this idea were based more on brute-force stochastic search rather than optimization theory, and the results were somewhat disappointing. Samuels' classic checkers-playing program [60] has been interpreted, in retrospect, as a kind of adaptive critic system. His adaptive “static position evaluator” served as a kind of Critic. More recently, Tesauro's master-class backgammon program [24,55] has clearly demonstrated that reinforcement learning systems can generate intelligent behavior.

[0079] In 1968, it was argued [20] that reinforcement learning could be used as a foundation for understanding intelligence in the brain. The concept of reinforcement learning (illustrated in FIG. 4) in the modern way was described. In this concept, a system learns a strategy of action that tries to maximize the long-term future expected value of utility, following the concepts of Von Neumann. For the first time, it was pointed out that a machine could be built to perform reinforcement learning simply by trying to approximate dynamic programming as formulated by Howard [28].

[0080] In 1972, the Werbos Harvard Ph.D. thesis proposal included a flowchart virtually identical to FIG. 2, with a detailed discussion of how to adapt the components of the system. However, the actual thesis focuses on a rigorous, generalized formulation of backpropagation as such, along with a demonstration of its effectiveness in system identification and political forecasting and a minimal discussion of neural networks [26,61]. (The thesis is reprinted in its entirety in [26].) Originally there were concerns about the learning speed of HDP, as one scales up to larger problems [8,22]. From 1977 to 1981, the ideas of HDP and of two more sophisticated methods—DHP and GDHP—were published addressing these scaling problems [25,62,63]. In particular. [63] discussed the idea of using a neural network trained by backpropagation in order to approximate the J function of dynamic programming.

[0081] In the meantime, on a totally independent basis, Bernie Widrow [23] published the first working, successful implementation of an adaptive critic, using a neural network as a Critic network. His 1973 paper is the original source of the term “Critic.” However, the Critic adaptation method that he used was not a true real-time learning method and has not been used elsewhere. In 1983, Barto, Sutton and Anderson [21] published a classic paper which became extremely famous in the 1990's. They implemented an adaptive critic system consisting of two adaptive elements or neurons. One element was a Critic, trained by a method which Sutton calls a Temporal Difference (TD) method. This method was developed on the basis of intuitive arguments about how to improve upon Widrow's algorithm. The method used in that paper was a very narrow special case of HDP. Later, in a famous 1988 paper [64], Sutton [64] expanded the concept of temporal difference methods to include “TD(λ)” for λ in the interval [0,1]. TD(1) essentially reproduces die old Widrow method, and TD(0) a generalization of the method in [21]. Very few papers actually use TD(λ) for λ other than zero (see 15,24,55). The other element of the system in [21] was a kind of Action network, trained by an associative reward/punishment algorithm called Arp, which is not a method for adapting Critics as such.

[0082] In January 1987 [65], a paper on ADP, GDHP, potential engineering applications and links to neuroscience was published, once again emphasizing the connection between reinforcement learning and dynamic programming. This paper led to major efforts by Barto and Sutton to follow up on this new connection to dynamic programming. Barto encouraged the early efforts in this area by Bertsekas [24], who, along with Tsitsiklis, has substantially enlarged the concept of TD. The theory behind TD his been expanded considerably by these and other efforts; however, the method itself remains a special case or subset (proper or otherwise) of HDP.

[0083] Offline simulation or “dreaming” has also been discussed as a way of improving the performance of ADP systems (as previously discussed in [65]). (Current evidence from psychology is consistent with this as a partial interpretation of dreams in human beings as well [66].) This was the basis of the “Dyna” architecture first presented by Sutton in [22]. This kind of offline learning is extremely important; however, the present invention is more suited to real-time systems.

[0084] The discussions of 1987 also contributed to the arrangement of an NSF workshop on neurocontrol, held in New Hampshire in 1988, which resulted in [22]. That workshop and book helped to bring together neurocontrol and adaptive critics as organized fields of research. FIG. 3.1 of that book illustrated the modem concept of reinforcement learning in a more picturesque way, equivalent to FIG. 4 above. Earlier forms of “reinforcement learning” that did not address the mathematical problem of maximizing utility (or minimizing cost) over time [67] did not help to carry through the control strategy given above with reference to equation (12).

[0085] More recently, “Action-Dependent Adaptive Critics” (ADAC) have been discussed. This generally includes Q-learning, ADHDP, ADDHP and ADGDHP [8,14], all of which are closely related. In fact, many of the “new” designs for “extended” or “modified” or “policy” Q-learning are actually implementations of ADHDP, which was reviewed at length in 1992 [14]. That book also reported a successful implementation at McDonnell-Douglas in the manufacturing of composite parts. It reported a successful simulated application to the control of damaged aircraft, which later led to a large related pro gram at NASA Ames. All of those methods use Critic networks which input not only R(t) but also u(t). Instead of approximating the J function of dynamic programming, they approximate a related function. When adapting a Critic where the controller or Action network is held fixed, the action-dependent methods for Critic adaptation reduce to their older, action-independent equivalents.

[0086] Another family of reinforcement learning methods called Alopex has also receievd attention. The most effective version of Alopex appeals to have been developed by Tzanakou of Rutgers. In that version, the Critic network is an action-dependent system, adapted by a procedure quite similar to ADHDP or Q-learning. Most of Klopf's final designs, in collaboration with Leemon Baird, have the same property. Although Q-learning and ADHDP appear less relevant, the two-sample method [50] discussed later is a significant exception.

[0087] In certain special cases, the same mathematics that underlie FIG. 2 can be implemented without a Model network as such. For example, in a plant with a single state variable, it is enough to know the sin of (∂R(t+1)/∂u(t)) in order to compute the derivatives of Ĵ(t+1) with respect to the weights, to within a scalar factor which effectively just modifies the learning rate α. This kind of effect has permitted the development of architectures similar to FIG. 2, but simpler, with comparable performance in certain special cases. (For example, see Lewis [68] and by Berenji [69].)

[0088] More recently, Error Critics and a multiple-model hierarchical decision system, both motivated by recent findings from neuroscience [8-10,13,32,46], have been discussed. The hierarchical decision system is a higher-level learning design, which essentially requires ordinary Critic networks as subsystems.

[0089] Issues of stability and convergence have been a major concern of the adaptive critic community for a long time. Nevertheless, formal results about MBAC systems have been relatively sparse. There have been the Ph.D. theses of Prokhorov [15] and of Landelius [43], and unpublished results from Richard Saeks of Accurate Automation. Sacks' work focuses on whole-system stability issues, rather than Critic adaptation as such. Whole-systems stability is the major concern in [15] and [43], but they also consider Critic adaptation as such. Landelius mainly proves stability and convergence for his own Critic implementations. Prokhorov [15] proves a whole-systems stability result for an offline learning version of HDPG+BAC in the deterministic case. He also reports considerable empirical testing of many approaches, as will be discussed below.

[0090] Tsitsiklis and Van Roy [73] have reviewed a large number of papers from the literature on lookup-table Critics, showing many specific examples where TD—a special case of HDP—becomes unstable. However, they prove that TD with lookup tables is stable, under certain restrictive assumptions. The key assumption is that the training examples x(t) must be taken from the probability distribution implied by the Markhov chain being studied. Stability against arbitrary training sequences is not guaranteed, even in the case of lookup table Critics.

[0091] In addition to these important but sparse theoretical results, considerable practical energy has gone into the effort to make these methods converge, for a variety of different test problems. This has been a very important practical issue, even in the context of offline learning. Some of this empirical work focused on the interplay between Action networks, Model networks and Critics. Narendra in [41] addresses related issues. But this paper focuses on Critic networks as such.

[0092] Among the tricks that are relevant to Critics as such, which have proven to be very useful in obtaining convergence, are:

[0093] (1) “Shaping,” [7,8,453]. In shaping, one first adapts a network to solve one class of control problems, and then uses the resulting weights as the initial values of the weights of a network trained to solve similar but more difficult problems. When this method is applied over and over again, it provides a way of doing “step-by-step learning” (not the same as real-time or incremental learning!) analogous to the way that humans really learn to perform difficult tasks.

[0094] (2) Interest rate management. In this approach, one starts out with very large values of r, and then one gradually decreases r to zero. (Usually r=0 represents our true values, if the utility function is properly articulated [12,74].) Strictly speaking, this is just a special case of shaping, because different values of r represent different optimization problems. Large values of r represent shorter-term optimization problems, which are usually easier to solve than the long-term problems.

[0095] (3) Utility function management. This can also be viewed as a form of shaping, in a way. In some cases, the choice of utility functions has even been used as way of cheating, of injecting prior knowledge about the desired controller into the adaptive system. However, for complex reasons [10,11,674], it does make sense to think of a hybrid man-machine control system, in which the human (as the upper controller) passes on something like his learned J function to the lower controller (the machine), which treats those inputs as its fundamental values (U). As a practical example—when using gradient-based MBAC methods to search for optimal energy-saving “chaos controllers,” it makes sense to use a relatively flat cost function U, representing energy consumption, in some large acceptable region of state space, and then to add a gradually rising penalty function (e.g., quadratic) for points in some “buffer zone” between that region and regions which are truly dangerous.

[0096] These three tricks have been successful in a very wide range of studies. Also, many researchers in ADP for discrete state space reported similar results, in discussions organized by Sridhar Mahadevan of the University of South Florida, in connection with an NSF workshop in April 1996 on reinforcement learning oriented towards A1.

[0097] It is an object of the present invention to provide a more stable control system design than previously available in the prior art.

[0098] This and other objects of the present invention are achieved by providing a new adaptive control system, implemented in either hardware or software (or as a hybrid of both). The control system (and its corresponding method) is designed to enhance the stability of the response of a controlled application (e.g., a factory, an airplane, a missile interception system. or a financial system).

[0099] A more complete appreciation of the invention and many of the attendant advantages thereof will be readily obtained as the same becomes better understood by reference to the following detailed description when considered in connection with the accompanying drawings, wherein:

[0100] FIG. 1 is a block diagram of an Indirect Adaptive Control System;

[0101] FIG. 2 is a block diagram of a combined HDP and BAC control system;

[0102] FIG. 3 is a hierarchy showing the early origins of the adaptive critic field; and

[0103] FIG. 4 is a block diagram of a reinforcement learning system.

[0104] A computer implements the method of the present invention. A computer housing houses a motherboard which contains a CPU, memory (e.g., DRAM, ROM, EPROM, EEPROM, SRAM, SDRAM, and Flash RAM), and other optional special purpose logic devices (e.g., ASICs) or configurable logic devices (e.g., GAL and reprogrammable FPGA). The computer also includes plural input devices, (e.g., a keyboard and mouse), and a display card for controlling monitor. In addition, the computer system further includes a floppy disk drive; other removable media devices (e.g., compact disc, tape, and removable magneto-optical media; and a hard disk, or other fixed, high density media drives, connected using an appropriate device bus (e.g., a SCSI bus, an Enhanced IDE bus, or a Ultra DMA bus). Also connected to the same device bus or another device bus, the computer may additionally include a compact disc reader, a compact disc reader/writer unit or a compact disc jukebox. The compact disc may be housed in a caddy or inserted directly into CD-ROM drives which do not require caddies. In addition, a printer also provides printed listings of control parameters (matrices and weights) used to implement the system.

[0105] As stated above, the system includes at least one computer readable medium. Examples of computer readable media are compact discs, hard disks, floppy disks, tape, magneto-optical disks, PROMs (EPROM, EEPROM, Flash EPROM), DRAM, SRAM, SDRAM, etc. Stored on any one or on a combination of computer readable media, the present invention includes software for controlling both the hardware of the computer and for enabling the computer to interact with a human user. Such software may include, but is not limited to, device drivers, operating systems and user applications, such as development tools. Such computer readable media further includes the computer program product of the present invention for implementing a stablized controller. The computer code devices of the present invention can be any interpreted or executable code mechanism, including but not limited to scripts, interpreters, dynamic link libraries, Java classes, Active X components, and complete executable programs.

[0106] Given this state of the art, many control engineers have questioned the practical value of adaptive control and real-time learning in general. However, extension of prior methods has practical engineering importance, for at least three reasons:

[0107] (1) Even when the stochastic optimization problem is solved offline, improved convergence can be of great practical value. A fully converged solution for J already implies a stability guarantee for the resulting feedback controller; however, the designers of such controllers would benefit from having tools which are more certain to provide these kinds of solution.

[0108] (2) Adaptation of parameters or networks “on the fly”—beyond the scope of substantive prior knowledge or models—is essentially the only tool available to address type-three uncertainty. Like biological brains., learning on the fly should be combinable together with application-specific feedback loops (which can also be tuned in real time), in order to get the best of both worlds.

[0109] (3) The extensions of adaptive control may improve stability and performance over time, and thereby enable this class of designs to become more competitive in some practical applications. These methods have inherent advantages over explicit multi-period optimization, such as reduced computational cost and correct treatment of noise in the nonlinear case. It is also possible to combine the two approaches, by using the ADP methods to estimate future costs or payoffs beyond the horizon of the n-step lookahead [15,33]. (Note that this kind of combination is quite different from the “rollout” methods discussed by Bertsekas [24], which use search trees rather than derivative-based methods like BTT.)

[0110] Thus, the present invention is directed to better approximate a true adaptive, control design that meets the ideal requirements for a universal controller as discussed above. The present invention is then extended to nonlinear, stochastic learning controllers that are applied to linear plants. The present invention draws upon the close parallels between a standard Indirect Adaptive Control (IAC) design, illustrated in FIG. 1, and the HDP+BAC design illustrated in FIG. 2. The design has three advantages. First, in linear quadratic optimization problems, it is well known that the solution to Bellman's equation is simply quadratic in the state of the system. Therefore, by choosing the “Critic network” to be a general quadratic function of the state, some set of weights will yield an exact solution of the Bellman equation. (This set must be verified empirically.)

[0111] Second, if a deterministic plant of this sort is controllable at all, then there will always exist some control law that sends the tracking error to zero, at an exponential rate. This guarantees that there exists a strategy or policy of action (control rule) for which the sum of future cost will be finite, even with r=U₀=0. Furthermore, dynamic programming always yields the strategy of action of minimum total future costs, when a minimum exists. That is, if a minimum does exist, the dynamic programming solution will do as well or better.

[0112] Third, Landelius [43] has already proven whole-system stability results for a particular implementation of HDP+BAC and for other methods to be discussed in greater detail below. Unfortunately, his new ways of implementing these methods—while very reasonable from a practical point of view—are technically not parallel, distributed systems. They are not true adaptive control designs in the spirit of Narendra, because they do not meet the tight constraints on computational cost and structure that are required for neural networks. Nevertheless, they are an important achievement. They go far beyond the earlier, similar work of Bradtke [44] in addressing total-system stability issues and indicate that this line of research can be successful.

[0113] In summary, if a parallel distributed version of HDP does in fact make the Critic converge to the dynamic programming solution, then some corresponding form of HDP+BAC is expected to provide a universal adaptive controller for the linear case. Such a structure would be a true adaptive control design, in the spirit of Narendra, because all of the concurrent adaptation laws are low-cost, single-pass adaptation laws. That is consistent with the local, parallel computing found in the brain—a style of computing which also makes it possible to build very high throughput dedicated chips.

[0114] Even for the fully observable linear deterministic case, it will take substantial effort to prove total system stability for this class of design, accounting for the coupled dynamics of all four parts of the overall system—the plant, the Action network, the Model network and the Critic. It took many years before Narendra and others [1-3] proved restrictive theorems about total system stability even for the SISO state, for systems with only three major components (i.e plant, Model. Action).

[0115] The first, most critical step in this process is to prove stability for the Critic adaptation as such, holding the rest as fixed, for a Critic adaptation rule that is also applicable to the nonlinear, stochastic case. However, none of the previously existing Critic adaptation rules were able to pass that test, for the definition of stability discussed below. The present invention provides new adaptation rules that do.

[0116] Critic adaptation methods are a way to derive Liapunov functions in a computer-based numerical fashion. The J function of dynamic programming is guaranteed to be a Liapunov function, at least for deterministic systems. This connection points towards another new direction for possible research [49]. One way of implementing this connection would be to use analytically-derived Liapunov functions as initial values for Critic networks. For example, a Liapunov function developed for rigid-body robotics could be tuned further, adaptively, using a method like HDP, in order to derive a more accurate Liapunov function for a somewhat flexible robot. By successive tuning, one could derive Critics suitable for the control of much more flexible robots.

[0117] In summary, HDP, DHP and GDHP are the three Critic adaptation methods relevant to the objects discussed herein. TD, as a subset (proper or improper) of HDP, is automatically included in this analysis. (The two-sample method is also discussed herein.) HDP and GDHP are methods for adapting either neural networks or other parametrized systems Ĵ(R(t), W)—with an array W of parameters or weights—so as to approximate the J function of dynamic programming. The J function of dynamic programming satisfies the extended Bellman equation (equation 12).

[0118] In FIG. 2, the Action network or controller implements a mathematical function u(R(t)). In the language of dynamic programming [24], this function is called a “policy” (π) or a strategy of action. The key theorem of dynamic programming is that the choice of actions u implied by equation 12 represents the best policy of action. In other words, for the case where U₀=0, it yields the policy π such that: $\begin{array}{cc}J\left(\underset{\_}{R}\left(t\right)\right)=\underset{\pi }{\mathrm{Max}}{J}^{\pi }\left(\underset{\_}{R}\left(t\right)\right),& \left(13\right)\end{array}$

[0119] where J^π is the expected value under policy π: $\begin{array}{cc}⟨\sum _{k=0}^{\infty }U\left(\underset{\_}{R}\left(t+k\right),\underset{\_}{u}\left(t+k\right)\right)/{\left(1+r\right)}^k⟩& \left(14\right)\end{array}$

[0120] For a fixed policy of action π, simple algebraic manipulations [24,44] yield an equation parallel to equation 12:

[0121] J^π((t))=U((t), ^π(t))+(1/(1+r))<J^π((t+1))>−U₀^π,  (15)

[0122] where U₀^π is the expected value of U across all future times under policy π. In finite-horizon problems with r=0, the best policy is the one which maximizes U₀^π. In the absence of “crossroads phenomena,” it can be obtained by choosing u(t) as in equation 12 for the maximal U₀^π. [28].

[0123] DHP attempts to adapt a network {umlaut over ({circumflex over (e)})}(R(t), W) so as to approximate the vector:

[0124] λ(t)=∇J((t)),  (16)

[0125] i.e., $\begin{array}{cc}{\lambda }_i\left(t\right)=\frac{\partial }{\partial R_i}J\left(\underset{\_}{R}\left(t\right)\right)& \left(17\right)\end{array}$

[0126] These λ_ivariables correspond to Lagrange multipliers [4,5] in deterministic optimal control, and to the “costate variables” in Pontryagin's approach. In economics, such Lagrange multipliers are associated with price signals, or, more precisely, with “shadow prices.” Thus, from an economist's viewpoint, DHP may be regarded as a distributed, adaptive system for calculating prices under conditions of uncertainty.

[0127] The recurrence relation for λ may be derived simply by differentiating equation 12. In the simplest, most direct representation, this yields [8, ch. 13]: $\begin{array}{cc}\begin{array}{c}{\lambda }_i\left(\underset{\_}{R}\left(t\right)\right)=\frac{\partial J\left(\underset{\_}{R}\left(t\right)\right)}{\partial R_i\left(t\right)}=\frac{\partial U(\underset{\_}{R}\left(t\right),\underset{\_}{u}\left(\underset{\_}{R}\left(t\right)\right)}{\partial R_i\left(t\right)}+\\ \frac{1}{1+r}⟨\frac{\partial J\left(\underset{\_}{R}\left(t+1\right)\right)}{\partial R_i\left(t\right)}⟩\\ =\frac{\partial U(\underset{\_}{R}\left(t\right),\underset{\_}{u}\left(\underset{\_}{R}\left(t\right)\right)}{\partial R_i\left(t\right)}+\sum _j\frac{\partial U\left(\underset{\_}{R},\underset{\_}{u}\right)}{\partial u_j}·\frac{\partial u_j\left(\underset{\_}{R}\left(t\right)\right)}{\partial R_i\left(t\right)}+\\ \frac{1}{1+r}\{\sum _j⟨\frac{\partial J\left(\underset{\_}{R}\left(t+1\right)\right)}{\partial R_i\left(t+11\right)}·\frac{\partial R_j\left(t+1\right))}{\partial R_i\left(t\right)}⟩+\\ \sum _{j,k}⟨\frac{\partial J\left(\underset{\_}{R}\left(t+1\right)\right)}{\partial R_i\left(t+11\right)}·\frac{\partial R_j\left(t+1\right))}{\partial u_k\left(t\right)}·\frac{\partial u_k\left(t\right)}{\partial R_j\left(t\right)}⟩\}\\ =\frac{\partial U\left(t\right)}{\partial R_i\left(t\right)}+\sum _j\frac{\partial U\left(t\right)}{\partial u_j}·\frac{\partial u_i\left(t\right)}{\partial R_i\left(t\right)}+\\ \frac{1}{1+r}\{\sum _j{\lambda }_j\left(t+1\right)(\frac{\partial R_j\left(t+1\right))}{\partial R_i\left(t\right)}+\\ \sum _k\frac{\partial R_j\left(t+1\right))}{\partial u_k\left(t\right)}·\frac{\partial u_k}{\partial R_j})\}\end{array}& \left(18\right)\end{array}$

[0128] Unfortunately, equation 18—unlike equation 12—no longer includes the recipe for choosing u. However, (as in equation 12) u should be chosen so as to maximize (or minimize) U(t)+<J(t+1)>/(1+r). For differentiable systems, this implies: $\begin{array}{cc}\begin{array}{c}0=\frac{\partial }{\partial u_i}\left(U\left(t\right)+\frac{1}{1+r}⟨J\left(t+1\right)⟩\right)\\ =\frac{\partial U\left(t\right)}{\partial u_i\left(t\right)}+\frac{1}{1+r}\sum _j⟨\frac{\partial J(\underset{\_}{R}\left(t+1\right)}{\partial R_j\left(t+1\right)}·\frac{\partial R_j\left(t+1\right)}{\partial u_i\left(t\right)}⟩\\ =\frac{\partial U\left(t\right)}{\partial u_i\left(t\right)}+\frac{1}{1+r}⟨{\lambda }_j\left(t+1\right)·\frac{\partial R_j\left(t+1\right)}{\partial u_i\left(t\right)}⟩\end{array}& \left(19\right)\end{array}$

[0129] Equation 18 is a stochastic version of the classical Pontryagin principle. In recent years, the mathematics community has shown interest in other, more formal stochastic versions of the Pontryagin principle [70]. For any particular policy π, the differentiation of equation 15 again yields equation 18, but without the additional requirement of equation 19.

[0130] As discussed herein, the Critic networks (Ĵ or {umlaut over ({circumflex over (e)})}) of HDP and DHP systems are adapted by gradient-based error minimization, as in classical backpropagation [26]. However, in the neural network community, many researchers prefer to use different types of adaptation, in order to perform basic tasks in supervised learning. (HDP previously has been detailed [8,17] in the general case so that the user can choose a preferred supervised learning system. DHP was similarly described in [8].) This is the reason why Landelius [43] was able to use an alternative supervised learning method, somewhat analogous to Recursive Least Squares (RLS), in his implementations of HDP, DHP, etc. (RLS is not a parallel distributed processing method, according to the strict standards of neural network research, but it comes close.) Unfortunately, in equation 10 of [8, ch. 13], there is a typographical error. The left-hand side should read “λ*(t)”.) The reader should also be warned that the pseudocode in [22] contains typographical errors. GDHP inherently requires a gradient-based learning approach.

[0131] In practice, the restriction to gradient-based learning is not as serious as it may seem at first, because many popular learning algorithms can be expressed in that form. For further examples and discussions of how to implement these methods in the general case, see the literature discussed above with reference to adapting Critics. Certain kinds of second-order information can also be used in true gradient-based parallel processing designs, as will be discussed later. Landelius has suggested that his ways of implementing DHP and ADDHP might provide an alternative route for achieving a low-cost universal adaptive control.

[0132] To adapt the Action networks, in the general case, one mainly uses the derivatives of U(t)+<J(t+1)> with respect to the weights in the Action network. One can calculate those derivatives by first calculating the derivatives with respect to u(t) (as in equation 19, implemented by backpropagation), and then using backpropagation again to propagate them through the Action network. For HDP and GDHP, see the discussion above of HDP+BAC, except that an extra term is added for the derivatives of U(R(t), u(t)) with respect to u(t). The present invention, however, focuses on the special case where U(t) is expressed simply as U(R(t)), so that this additional, term drops out. Even though equation 19 calls for the use of an expectation value, the same change of weights can be achieved, on average, by simply adapting the weights in response to the currently observed derivatives, as in the usual LMS learning and statistical sampling.

[0133] The present invention focuses on assessing the stability and accuracy of Critic adaptation in situations where the Action network or policy is held fixed, where the Model network is also held fixed, where U(t) does not depend on u(t), and where Critic adaptation is based on the traditional gradient-based versions of the Critic adaptation methods. In addition, the present invention assumes a completely observable plant, where R(t) and X(t) are the same. For simplicity, it is assumed that U₀=0 or that U₀is known apriori.

[0134] In this special case, there really is no Action network, in effect. Thus, the augmented Model network is defined ,is the combination of the Model network proper and the Action network, treated as one large network. The augmented Model network is written as:

[0135] (t+1)=((t), (t)),  (20)

[0136] where w(t) refers to a vector of random variables, as in [39]. (Statisticians more often use the notation e(t) or u(t) to refer to a vector of random numbers; however, this notation has led to severe confusion in parts of the control theory community, especially involving the definition of “ARMA” processes, which are explained at length in [26,29].)

[0137] In real-time adaptation in HDP, the Critic is adapted as follows at each time t:

[0138] (1) First obtain the vector R(t), simply by observing (or estimating) the current state of the plant.

[0139] (2) Next obtain a sample of R(t+1), either by simulating equation 20 or by waiting until time t+1 and obtaining a new observation.

[0140] (3) Update the weights W_ijby: $\begin{array}{cc}{W}_{\mathrm{ij}}\left(t+1\right)=W_{\mathrm{ij}}\left(t\right)-\alpha \frac{\partial E}{\partial W_{\mathrm{ij}}}& \left(21\right)\end{array}$

[0141] where α is a learning rate (as in equations 8 and 11), and

[0142] E=½(e(t))²,  (22a)

[0143] where:

[0144] e(t)=Ĵ((t+1))/(1+r)+U(t)−Ĵ((t))−U₀  (22b)

[0145] This is not a complete description of the problem, because it does not specify how to compute the derivatives that appear in equation 21. This is not just a matter of computation; there is also a need to specify which partial derivatives or gradient are used. In HDP proper, the system calculates: $\begin{array}{cc}\frac{\partial E}{\partial W_{\mathrm{ij}}}=2e\left(\frac{\partial E}{\partial W_{\mathrm{ij}}}\right)=2e\left(-\frac{\partial }{\partial W_{\mathrm{ij}}}\hat{J}\left(\underset{\_}{R}\left(t\right),W\right)\right)& \left(23\right)\end{array}$

[0146] This calculation treats J(R(t)) in equation 22 as a function of W, but does not account for the effect of W on J(R(t+1))! If equation 12 were solved using an extension of the Galerkin method [27], the weights would be adapted by use of the complete total gradient: $\begin{array}{cc}\frac{\partial E}{\partial W_{\mathrm{ij}}}=2e\left(-\frac{\partial }{\partial W_{\mathrm{ij}}}\hat{J}\left(\underset{\_}{R}\left(t\right),W\right)+\frac{1}{1+r}\frac{\partial }{\partial W_{\mathrm{ij}}}\hat{J}\left(\underset{\_}{R}\left(t+1\right),W\right)\right)& \left(24\right)\end{array}$

[0147] Combining equations 21, 22 and 23, pure HDP may be written as:

[0148] ΔW=−2αe(∇_pe)=2αe∇Ĵ((t))  (25)

[0149] and the Galerkinized version of HDP, HDPG, may be written as:

[0150] ΔW=−2αe(∇_te)=2αe(∇Ĵ((t))−(1/(1+r))∇Ĵ((t+1))),  (26)

[0151] where ∇_prefers to the partial gradient with respect to the weights, and where ∇_trefers to the total gradient. All of these derivatives can be obtained very efficiently by use of backpropagation, using local parallel calculations and so on [8,26].

[0152] In DHP, essentially the same three steps are followed, but the error, E, is defined as: $\begin{array}{cc}E={\left(\underset{\_}{e}\left(t\right)\right)}^2=\sum _i{\left({\lambda }_i*\left(t\right)-{\hat{\lambda }}_i\left(t\right)\right)}^2,& \left(27\right)\end{array}$

[0153] where the target value λ_i*(t) is calculated by: $\begin{array}{cc}{\lambda }_i^{*}\left(t\right)=\frac{\partial U\left(t\right)}{\partial R_i}+\frac{1}{1+r}\sum _j{\hat{\lambda }}_j\left(t+1\right)·\frac{\partial f_j\left(\underset{\_}{R(}t\right),\underset{\_}{w}\left(t\right))}{\partial R_i\left(t\right)}& \left(28\right)\end{array}$

[0154] Note that {circumflex over (λ)}_i(t) is simply the i-th component of the vector λ(R(t),W) calculated by the Critic network, and that {circumflex over (λ)}_j(t+1) is the j-th component of {umlaut over ({circumflex over (e)})}(R(t+1),W). Equation 28 corresponds to equation 10 of [8,ch.13] in the case where there is no Action network (and where the typo on the left is fixed). However [8,ch. 13] explicitly indicates how the derivatives may be calculated efficiently, using “dual subroutines” or “dual functions” in order to implement backpropagation. Note as well that the backpropagation equations for equation 28 require that the vector w(t) be known. When R(t+1) was obtained from simulation, w(t) was already known, but when R(t+1) is observed directly, w is deduced efficiently only if the Model network has a special kind of structure. (For example, the Model network could be based on the usual assumption that R_i(t+1)={circumflex over (R)}_i(t+1)+w_i, or it could be based on a structure like the Stochastic Encoder/Decoder, Predictor of [8,ch.13].)

[0155] Here, as with HDP, either a pure version of DHP, based on ∇_pE, or a Galerkinized version, based on ∇_tE, may be chosen. In pure DHP: $\begin{array}{cc}\Delta W=-2\alpha \sum _k{e}_k{\nabla }_p{e}_k=2\alpha \sum _k{e}_k\nabla {\hat{\lambda }}_k\left(\underset{\_}{R}\left(t\right)\right)=-2\alpha \left({\underset{\_}{e}}^T{\nabla }_p\underset{\_}{e}\right)=2\alpha \left({\underset{\_}{e}}^T\nabla \underset{\_}{\stackrel{\widehat{¨}}{e}}\left(t\right)\right)& \left(29\right)\end{array}$

[0156] In DHPG: $\begin{array}{cc}\begin{array}{c}\Delta W_{\mathrm{ij}}=-2\alpha \sum _k{e}_k\left(\frac{\partial {\lambda }_k^{*}\left(t\right)}{\partial W_{\mathrm{ij}}}-\frac{\partial \widehat{{\lambda }_k}\left(t\right)}{\partial W_{\mathrm{ij}}}\right)\\ =-2\alpha \sum _k{e}_k(-\frac{\partial \widehat{{\lambda }_k}\left(t\right)}{\partial W_{\mathrm{ij}}}+\\ \frac{1}{1+r}\sum _v\frac{\partial f_v\left(R\left(t\right),\underset{\_}{w}\left(t\right)\right)}{\partial R_k}·\frac{\partial \widehat{{\lambda }_v}\left(t+1\right)}{\partial W_{\mathrm{ij}}})\\ =2\alpha \sum _k{e}_k\frac{\partial \widehat{{\lambda }_k}\left(t\right)}{\partial W_{\mathrm{ij}}}-\frac{2\alpha }{1+r}\sum _{k,v}{e}_k\frac{\partial f_v}{\partial R_k}·\frac{\partial \widehat{{\lambda }_v}\left(t+1\right)}{\partial W_{\mathrm{ij}}}\end{array}& \left(30\right)\end{array}$

[0157] This can be represented as: $\begin{array}{cc}\Delta W=-2\alpha \sum _k{e}_k{\nabla }_t{e}_k=2\alpha \left({\underset{\_}{e}}^T\nabla \underset{\_}{\stackrel{\widehat{¨}}{e}}\left(t\right)\right)-\frac{2\alpha }{1+r}{\underset{\_}{e}}^T{F}^T\left(t\right)\nabla \underset{\_}{\stackrel{\widehat{¨}}{e}}\left(t+1\right),& \left(31\right)\end{array}$

[0158] where F(t) is the Jacobian with respect to R(t) of the function f(R(t),w(t)), and where all of the gradients are gradients with respect to the array of weights W.

[0159] Finally, in GDHP, the same three steps are used as with HDP in order to adapt the weights in the Critic network J(R(t),W), but equation 22a is replaced with the error function:

[0160] E=Ω₀(e₀²(t))+^TΩ(t),  (32)

[0161] where Ω₀is some nonnegative scalar and Ω is some positive definite matrix, where e₀is defined to be the same as “e” as defined in equation 22b, and where e(t) is essentially the same vector as in equation 27. (The definition for e(t) changes because a {umlaut over ({circumflex over (e)})} output available is not available when using a critic which outputs Ĵ. Thus to define e(t) here, still using equations 27 and 28, each reference to {umlaut over ({circumflex over (e)})} is replaced by a reference to ∇Ĵ.) This is already sufficient to define GDHP, in a purely mathematical sense. Again, of course, there are two versions of the method, one a pure version and the other a Galerkinized version, depending on whether one adapts the weights in proportion to ∇_pE or ∇_tE.

[0162] GDHP is far more difficult to implement than HDP or DHP. Therefore, most researchers in this area focus their efforts on studying HDP and DHP, which already present major challenges and opportunities. Only Wunsch and Prokhorov have claimed a working implementation of GDHP[14]. To calculate the derivatives of E with respect to the weights, directly, second-order backpropagation is used; the theory of how to do this is given in [8,ch. 10], and pseudocode may be found in [71]. Indirect approaches may or may not converge to the same equilibrium. Even though GDHP is more difficult to implement than HDP or DHP, it is still a local, parallel neural network design, O(N) in cost in the same way as the other.

[0163] In [65], it was proposed that 92 should be a diagonal matrix, in effect. From a neural network viewpoint, that is a reasonable restriction. Herein, only two special cases are considered: (1) Ω=1; (2) Ω as any positive symmetric real matrix. In both cases, Ω₀=0. Strictly speaking, then, the conclusions below about GDHP apply only to those special cases of GDHP. But the general form of GDHP is a hybrid, in effect, of HDP (which minimizes e₀²) and of those special cases. It is difficult to believe that any form of GDHP could be unconditionally stable if those special cases and HDP both fail to achieve unconditional stability.

[0164] Using HDP, DHP or GDHP as methods for continuous-time adaptation appears to unique. There are several reasons for this. First, in real-world computing, digital or analog, there always are time delays, especially when complex nonlinear calculations are required. Continuous-time equations that assume that special-purpose analog chips are used and that no such delays exist would not accurately represent the actual behavior of a real-world computational system. Second, there are fundamental theoretical reasons why higher-order intelligent systems must contain discrete “clocks” [32].

[0165] Third, there is evidence from a host of respected neuroscientists (e.g. Llinas and Barry Richmond) that indicates that there are surprisingly precise “clocks” all over the brain. For example, there are general, nonspecific synchronized signals which enter the cerebral cortex, and create alternating “windows” of 30-50 ms in which the dynamics of the entire circuit appear to behave differently. Fourth, a discrete-time representation tends to simplify the communication and presentation of complex designs. Fifth, one can normally approximate the behavior of continuous-time systems with a series of discrete-time systems, but not vice-versa. Of course, the continuous time version of the Bellman equation—sometimes called the Hamilton-Jacobi-Bellman equation—has been well know for decades [4].

[0166] On the other hand, Narendra [1] and Grossberg [72] and others have proved a number of very important stability theorems, using a combination of Liapunov stability analysis and continuous-time representations. There is a connection between Liapunov theory and continuous-time representation.

[0167] For HDP, the discrete-time case and the continuous-time case do not appear to be very different. (Those cases are not likely to be very different for DHP or GDHP, either.)

[0168] The continuous-time version of equation 15 is:

[0169] ∂_tJ=−U+U₀+rJ  (33)

[0170] The continuous-time version of equation 21 is: $\begin{array}{cc}{\partial }_t{W}_{\mathrm{ij}}=-\alpha \frac{\partial }{\partial W_{\mathrm{ij}}}E=-2\alpha e\frac{\partial }{\partial W_{\mathrm{ij}}}e,& \left(34\right)\end{array}$

[0171] where:

[0172] e=∂_tĴ+U−U₀−rJ  (35)

[0173] For the continuous time version of pure HDP, HDPC, the partial gradient is given by: $\begin{array}{cc}\frac{\partial e}{\partial W_{\mathrm{ij}}}=-\frac{\partial \hat{J}}{\partial W_{\mathrm{ij}}}& \left(36\right)\end{array}$

[0174] (Multiply this by (1−r) has no real effect on equation (34) since it simply multiplies the arbitrary fixed scalar α by another fixed positive scalar.) For the continuous-time version of HDPG, HDPGC, the full gradient is given by: $\begin{array}{cc}\frac{\partial e}{\partial W_{\mathrm{ij}}}=\frac{d}{dt}\left(\frac{\partial \hat{J}}{\partial W_{\mathrm{ij}}}\right)-r\frac{\partial \hat{J}}{\partial W_{\mathrm{ij}}}& \left(37\right)\end{array}$

[0175] where the total derivative with respect to time is computed as if W were constant. In other words, the system defines: $\begin{array}{cc}\frac{d}{dt}\left(\frac{\partial \hat{J}}{\partial W_{\mathrm{ij}}}\right)=\sum _k\left(\frac{{\partial }^2\hat{J}}{\partial W_{\mathrm{ij}}\partial R_k}{\stackrel{.}{R}}_k\right)& \left(38\right)\end{array}$

[0176] which could be calculated in a number of ways (such as direct observation of the time-changes of the network under fixed weights or the use of second-order backpropagation).

[0177] The continuous-time version of DHP is based on the gradient with respect to R of equation 33:

[0178] ∂_tλ=−∇U+r  (39)

[0179] The continuous-time version of equation 21, for E defined as in equation 27, is: $\begin{array}{cc}{\partial }_t{W}_{\mathrm{ij}}=-\alpha \frac{\partial }{\partial W_{\mathrm{ij}}}\left|{\underset{\_}{e}}^2\right|=-2\alpha \sum _k{e}_k\frac{\partial }{\partial W_{\mathrm{ij}}}{e}_k& \left(40\right)\end{array}$

[0180] where: $\begin{array}{cc}{e}_k={\partial }_t\widehat{{\lambda }_k}+\frac{\partial U}{\partial R_k}-r\widehat{{\lambda }_k}& \left(41\right)\end{array}$

[0181] For the continuous-time version of pure DHP, DHPC, equation 40 is completed by defining: $\begin{array}{cc}\frac{\partial }{\partial W_{\mathrm{ij}}}{e}_k=-\frac{\partial }{\partial W_{\mathrm{ij}}}\widehat{{\lambda }_k}& \left(42\right)\end{array}$

[0182] For DHPGC, the system defines: $\begin{array}{cc}\frac{\partial }{\partial W_{\mathrm{ij}}}{e}_k=\frac{d}{dt}\left(\frac{\partial \widehat{{\lambda }_k}}{\partial W_{\mathrm{ij}}}\right)-r\frac{\partial \widehat{{\lambda }_k}}{\partial W_{\mathrm{ij}}}& \left(43\right)\end{array}$

[0183] where the total derivative with respect to time in equation 43 is defined exactly in parallel with the similar equations 37 and 38.

[0184] Finally, in GDHPC and GHDPGC, the weights are adapted according to: $\begin{array}{cc}{\partial }_t{W}_{\mathrm{ij}}=-\alpha \frac{\partial }{\partial W_{\mathrm{ij}}}E,& \left(44\right)\end{array}$

[0185] where E is defined as in equation 32, and where e₀,e(t) and their derivatives are based on e, e(t) and their derivatives for HDPC and DHPC—exactly in parallel to the way that HDP and DHP are used in GDHP in the discrete-time case.

[0186] As discussed above, it is assumed that the Action network is held fixed, which results in an augmented model R(t+1)=f(R(t), w(t)). In working with the linear observable case, notation similar to that of Narendra and Annaswamy [1] is used. Thus the state vector will be denoted simply as x. Under these conditions, a linear deterministic MIMO plant obeys the following very simple dynamics:

[0187] ∂_t(t)=A  (45)

[0188] for the continuous-time case, or:

[0189] (t+1)=A(t)  (46)

[0190] for the discrete-time case. The linear stochastic MIMO plant obeys:

[0191] ∂_t=A  (47)

[0192] in the continuous-time case, where w is a continuous stochastic process [4,5]. In the discrete-time case, it obeys:

[0193] (t+1)=A(t)+(t),  (48)

[0194] where w is a random vector based on the multivariate normal distribution N(0,Q), which implies that:

[0195] <>=0  (49)

[0196] <^T>=Q  (50)

[0197] By using a minimization of tracking error as an immediate goal, which is a quadratic in the state vector, the cost function to be minimized is assumed to be quadratic (technically, a quadratic disutility function) according to:

[0198] U()=^TV  (51)

[0199] where r is assumed to be zero.

[0200] For many linear-quadratic optimization problems, it is well known that the true J function of dynamic programming is quadratic in the state vector x. This was crucial to the research strategy discussed above. Thus for all of the Ĵ Critic networks applied to these problems:

[0201] Ĵ()=^TC  (52)

[0202] For the {umlaut over ({circumflex over (e)})} Critics, the gradient of Ĵ with respect to x is given by:

[0203] ()=2C  (53)

[0204] In order to analyze stability properties, a weight error matrix M is given by:

[0205] M=C−C*  (54)

[0206] where C* is the correct Critic weight matrix, the matrix which satisfies the Bellman equation for these systems:

[0207] J((t))=U((t))+<J((t+1)>−U₀  (55)

[0208] which, upon substitution from 51 and 53, reduces to:

[0209] C*=V+A^TC*A  (56)

[0210] both in the deterministic case (where U₀=0) and in the stochastic case (where U₀=<w^TC*w>=Tr(C*Q), a scalar independent of x).

[0211] For the stochastic version of these problems, with r=0, the U₀term is not generally equal to zero. This affects HDP and its variants, but not DHP or GDHP (as implemented here). For simplicity, U₀has not been adapted here, because it does not affect the issues presented herein. In practical computing work where r=0, the usual practice is to start with r>0 and gradually send r towards zero, and see what the Action network converges to; the adaptation of U₀has not received serious attention. This is a point worth returning to, as the basic understanding becomes more rigorous.

[0212] In describing the effects of adaptation in the discrete-time case, infinitely small learning, rates have implicitly been assumed. Many techniques exist to reduce learning rates, in practice, when they are so large that they lead to instabilities. Those issues are certainly worthy of study. Here, however, the more basic question is whether these methods could be unconditionally stable, even if there were no danger of learning rates being set too high. This approach simplifies the analysis very greatly, because the adaptation of C can be represented as a continuous process represented as a differential equation. This, in turn, allows the use of the usual continuous-time Liapunov stability methods.

[0213] To derive the adaptation rule for C with HDPC and HDPGC in the linear case with r=0, equations 51 and 52 are substituted into equation 35, which yields:

[0214] e=∂_t(^TC)+^TV₀=^TC^TC^TV₀,  (57)

[0215] where “∂_t” is a standard notation from physics for differentiation with respect to time. For the deterministic case, equation 45 is substituted into equation 57, and since U₀=0 here, e is given by:

[0216] e=^TA^TC^TCA^TV  (58)

[0217] For HDPC, equation 52 is substituted into equation 36, remembering that the weights “W_ij” for the Critic shown in equation 52 are simply the matrix elements C_ij. This yields: $\begin{array}{cc}\frac{\partial e}{\partial C_{\mathrm{ij}}}=-x_i{x}_j& \left(59\right)\end{array}$

[0218] Finally, substituting equations 59 and 58 into equation 34, the adaptation rule for HDPC in the linear deterministic case becomes:

[0219] ∂_tC_ij=α(^TA^TC^TCA^TV)x_ix_j  (60a)

[0220] which can be expressed more compactly as:

[0221] ∂_tC=αe^T  (60b)

[0222] For HDPGC, it is easiest to compute the total gradient of e with respect to C by differentiating equation 58 directly:

[0223] ∇_te=^TA^T+A^T  (61)

[0224] Substituting equations 61 and 58 back into equation 38:

[0225] ∂_tC=−αe(^TA^T+A^T)  (62)

[0226] Equation 54 is then solved to get:

[0227] C=M+C*  (63)

[0228] Substituting this into equation 58, and subtracting out the terms containing C* (which cancel out due to the Bellman equation, as expressed in equations 55 and 56), e is rewritten as:

[0229] e=^TA^TM^TMA  (64)

[0230] Since the matrix C* in equation 63 is constant (i.e. unchanged by our adaptation process), substitution into equation 60 yields equations (65) and (66) for HDPC and HDPGC, respectively.

[0231] ∂_tM=αe^T  (65)

[0232] ∂_tM=−αe(^TA^T+A^T).  (66)

[0233] For all the varieties of HDP, the matrices C*, C and M are assumed to be symmetric, without loss of generality, because: (1) only the symmetric part of C or C* affects the function J(x); (2) the update rules for C and M all involve the addition of symmetric matrices only (which can be seen by simply looking at the right-hand-side of these update equations).

[0234] Continuing with the analysis for HDP and HDPG, substituting equations 51 and 52 into 22b, again with r=0, produces:

[0235] e=(t+1)^TC(t+1)+^TV^TC₀  (67)

[0236] where “x” refers to x(t). Substituting equation 46 into equation 67, and exploiting the fact that U₀=0 in the deterministic case, produces:

[0237] e=^TA^TCA^TV^TC^T(A^TCA+V−C)  (68)

[0238] For HDP, substituting equation 52 into equation 23, produces: $\begin{array}{cc}\frac{\partial E}{\partial W_{\mathrm{ij}}}=-2{\mathrm{ex}}_i{x}_j& \left(69\right)\end{array}$

[0239] Assuming infinitely small learning rates in equation 21, such that equation 21 is converted into a differential equation, the substituting from equation 69 into that equation yields: $\begin{array}{cc}{\partial }_t{C}_{\mathrm{ij}}=-\frac{1}{2}\frac{\partial E}{\partial W_{\mathrm{ij}}}={\mathrm{ex}}_i{x}_j& \left(70\right)\end{array}$

[0240] Although an apparent learning rate of α=½ was chosen (for simplicity's sake), in actuality, the choice of apparent leaning rate has no effect on the meaning of this equation. It only rescales the arbitrary time scale of this equation.

[0241] For HDPG, equation 52 is substituted into equation 24, and then the result is substituted into equation 21 in the same maimer, which results in:

[0242] ∂_tC=−e(A^TA^T−^T)  (71)

[0243] Substituting equation 63 into equation 68, and again exploiting the Bellman equation, produces:

[0244] =^T(A^TMA−M)  (72)

[0245] Substituting equation 63 into equation 70, produces equations (73) and (74) for HDP and HDPG, respectively, given by:

[0246] ∂_tM=e^T  (73)

[0247] ∂_tM=e(^T−A^TA^T).  (74)

[0248] For the linear stochastic case, equation 48 is substituted into equation 57, and recall that U₀=Tr(CQ) here, which yields:

[0249] e=(A)^TC(A)+^TV^TC(C*Q)=^TA^TCA^TA^TC^TCA^TC^TV^TC(C*Q)  (75)

[0250] For HDP, equation 52 is substituted into equation 23, as in the deterministic case. The adaptation rule for C is still equation 70, but the meaning of equation 70 is changed by the use of a new definition for e, the definition given in equation 75. However, taking statistical expectations, with respect to the random vector w, produces (by in effect exploiting equations 49 and 50):

[0251] <∂_tC>=<e^T>=^T<e>=^T(^TA^TCA^TV^TC^TA^TC^TCA^TC(C*Q)>)=^T(^TA^TCA^TV^TC((C−C*)Q))  (76)

[0252] For HDPG in the stochastic case, equation 75 is differentiated with respect to C:

[0253] ∇_te=A^TA^T+^TA^T+A^T+^T−^T  (77)

[0254] Substituting into equation 24 and then 21, produces:

[0255] ∂_tC=−e(A^TA^T+^TA^T+A^T+^T−^T)  (78)

[0256] The expectation value of equation 78 is somewhat nasty to calculate, because it involves products up to the fourth power in the random vector w. However, as is apparent from the definition of e in equation 75, four of the terns do not depend upon w; pulling out their contribution, and then exploiting the fact that the expectation value of odd powers of w will always equal zero, equation 78 is expanded to:

[0257] <∂_tC>=−(^T(A^TCA+V−C)−Tr(C*Q))(<A^TA^T+^TA^T+A^T+^T−^T>)−<(^TA^TC^TCA)(^TA^T+A^T)>−<(^TC)(A^TA^T+^T−^T)>  (79)

[0258] The middle line of this equation is somewhat tricky to expand. As an example, one of its terms, <(x^TA^TCw)wx^TA^T>, can be expressed in an Einstein-like convention, in which doubled subscripts are implicitly summed over, yielding:

[0259] <(x_iA^T_ijC_jkw_k)w_μx_λA^T_λv>=x_iA^T_ijC_jkQ_kμx_λA^T_λv  (80)

[0260] which represents the matrix:

[0261] QCA^TA^T  (81)

[0262] Using this kind of approach, equation 79 is further expanded to produce:

[0263] <∂_tC>=−(^T(A^TCA+V−C)−Tr(C*Q))(A^TA^T+Q−^T)−2QCA^TA^T−2A^TA^TCQ−Tr(CQ)(A^TA^T−^T)−<(^TC)(^T)>=−(^T(A^TCA+V−C)+Tr((C−C*)Q))(A^TA^T+Q−^T)−2QCA^TA^T−2A^TA^TCQ+Tr(CQ)Q−<(^TC)(^T)>  (82)

[0264] Finally, for the analysis of error dynamics, equation 63 is again substituted into equations 76 and 82, and the Bellman equation is exploited. For HDP, this produces:

[0265] <∂_tM>=^T(^T(A^TMA−M)+Tr(MQ))  (83)

[0266] For HDPG, this produces:

[0267] <∂_tM>=−(^T(A^TMA−M)+Tr(MQ))(A^TA^T+Q−^T)−2Q(M+C*)A^TA^T−2A^TA^T(M+C*)Q+Tr((M+C*)Q)Q−<(^T(M+C*))(^T)>  (84)

[0268] Again, without loss of generality, it is assumed that C, C* and M are all symmetric matrices. In order to derive the details of DHP and DHPG in the stochastic linear case, with r=0, equations 51 and 48 are substituted into equation 28, to deduce:

[0269] λ*=2V^T(t+1)  (85)

[0270] Substituting 53 and 48 into 85, produces:

[0271] λ*=2V^T(2C)((A)=2V+2A^TCA+2A^TC  (86)

[0272] Substituting 53 and 86 into 27, produces:

[0273] =λ*−=2(V^TCA^TC)  (87)

[0274] To derive the details for DHP in this application, the derivatives that make up tile partial gradient, are calculated simply by differentiating equation 87 directly: $\begin{array}{cc}\frac{\partial }{\partial C_{\mathrm{ij}}}{e}_k=2\left(-{\delta }_{\mathrm{ik}}{x}_j\right)!& \left(88\right)\end{array}$

[0275] From this, it holds that: $\begin{array}{cc}\left({\underset{\_}{e}}^T\frac{\partial }{\partial W_{\mathrm{ij}}}\underset{\_}{e}\right)=\sum _k{e}_k\frac{\partial }{\partial C_{\mathrm{ij}}}{e}_k=-2e_i{x}_j& \left(89\right)\end{array}$

[0276] Substituting equations 87 and 89 into equation 29, and converting equation 29 into a differential equation produces:

[0277] ∂_tC=−½(^T∇_p)=(V^TCA^TC)^T=(V+A^TCA−C)xx^T+A^TC^T  (90)

[0278] The expected value over the random vector w is easily calculated:

[0279] <∂_tC>=(V+A^TCA−C)^T  (91)

[0280] Note that the deterministic linear case is simply the special case where w=0; in that case, ∂_tC=<∂_tC> is still given by equation 91.

[0281] In order to derive the details for DHPG, the values of two of the terms which appear in equation 30 must be determined for the linear stochastic case: $\begin{array}{cc}\frac{\partial f_v}{\partial R_k}=A_{\mathrm{vk}}& \left(92\right)\\ \frac{\partial \widehat{{\lambda }_v}\left(t+1\right)}{\partial W_{\mathrm{ij}}}=\frac{\partial }{\partial C_{\mathrm{ij}}}{\left(2C\left(A\underset{\_}{x}+\underset{\_}{w}\right)\right)}_v=2{{\delta }_{\mathrm{iv}}\left(A\underset{\_}{x}+\underset{\_}{w}\right)}_j& \left(93\right)\end{array}$

[0282] Substituting into the more complicated summation in equation 30 produces: $\begin{array}{cc}\sum _{k,v}{e}_k\frac{\partial f_v}{\partial R_k}·\frac{\partial \widehat{{\lambda }_v}\left(t+1\right)}{\partial W_{\mathrm{ij}}}=2\sum _{k,v}{e}_k{A}_{\mathrm{vk}}{{\delta }_{\mathrm{ij}}\left(A\underset{\_}{x}+\underset{\_}{w}\right)}_j=2{\left(A\underset{\_}{e}\right)}_i{\left(A\underset{\_}{x}+\underset{\_}{w}\right)}_j& \left(94\right)\end{array}$

[0283] Substituting back into equation 30, and translating these update equations into a differential equation as before, produces: $\begin{array}{cc}\begin{array}{c}{\partial }_{i}C={\underset{\_}{\mathrm{ex}}}^T-A\underset{\_}{e}\left({\underset{\_}{x}}^T{A}^T+{\underset{\_}{w}}^T\right)\\ =\left(V+A^T\mathrm{CA}-C\right){\underset{\_}{\mathrm{xx}}}^T+A^{T}C{\underset{\_}{\mathrm{wx}}}^T-\\ A\left(V+A^T\mathrm{CA}-C\right){\underset{\_}{\mathrm{xx}}}^T{A}^T-{\mathrm{AA}}^{T}C{\underset{\_}{\mathrm{wv}}}^T{A}^T-\\ A\left(V+A^T\mathrm{CA}-C\right){\underset{\_}{\mathrm{xw}}}^T-{\mathrm{AA}}^{T}C{\underset{\_}{\mathrm{ww}}}^T\end{array}& \left(95\right)\end{array}$

[0284] This implies:

[0285] <∂_tC>=(V+A^TCA−C)^T−A(V+A^TCA−C)^TA^T−AA^TCQ  (96)

[0286] In order to analyze the error dynamics of pure DHP, equation 63 is substituted into equation 91, and the Bellman equation is used to derive:

[0287] <∂_tM>=(A^TMA−M)^T  (97)

[0288] For DHPG, the same kind of substitution into equation 96 yields:

[0289] <∂_tM>=(A^TMA−M)^T−A(A^TMA−M)^TA^T−AA^T(C+M)Q  (98)

[0290] Note that the update rules shown in equations 91 and 96-98 do not in general add a symmetric matrix to C or M! With pure DHP or with DHPG, it cannot safely be assumed that C or M is a symmetric matrix: even if the initial estimate of C* is symmetric, the updated estimates after learning will usually not be symmetric. This is a valid consequence of the original nonlinear Critic adaptation rule that these special cases are being derived from. It is a manifestation of a general problem in the nonlinear case: the problem that DHP does not force λ(R(t)) to be an integrable vector field, even though the gradient of any scalar field (like J(R(t)) should be integrable. This problem was one of the main motivations behind the development of GDHP, which guarantees integrability even in the nonlinear case.

[0291] GDHP is more complicated than HDP or DHP in the general nonlinear case; however, in the linear stochastic case, with Ω=₀=0 (which allows us to ignore U₀), it reduces to something relatively simple. C might or might not be symmetric, for the sake of generality. For the case where Ω₀=0 and Ω=1, equation 32 becomes exactly the same as equation 27, except that alternative sources are used for the vectors λ* and {umlaut over ({circumflex over (e)})}. These vectors are defined by using ∇J in place of equation 53. Differentiating J as given in equation 52, produces:

[0292] λ()=(C+C^T)  (99)

[0293] Aside from this replacement, the farther calculations for GDHP are identical to those for DHP. When derivatives of E or e are calculated with respect to the weights C_ij, equation 99 simply replaces the previous (DHP) value of (∂E/∂C_ij) by the average ½(∂E/∂C_ij)+(∂E/∂C_ji)), because each weight C_ijappears both in C (as before) and in C^T. The net effect of this is that all of the adaptation equations for C, and dynamic equations for M, are exactly the same as for DHP and DHPG, except that every update is divided by two, and added it to its transpose. (Also, the division by two can be skipped, since it simply multiplies everything by a scalar factor that can be absorbed into the learning rate.)

[0294] Thus for GDHP in the stochastic case, equation 91 is replaced by:

[0295] <∂_tC>=(V+A^TCA−C)^T+^T(V+A^TCA−C)  (100)

[0296] Its error dynamics are:

[0297] <∂_tM>=(A^TMA−M)^T+^T(A^TMA−M)  (101)

[0298] For GDHPG, the error dynamics follow the symmetrization of equation 98:

[0299] <∂_tM>=(A^TMA−M)^T+^T(A^TMA−M)−A(A^TMA−M)^TA^T−AA^T(C*+M)Q−A^T(A^TMA−M)A^T−Q(C*+M)AA^T  (102)

[0300] Strictly speaking, the asymmetry problem with DHP could have been patched up by using a neural network trick called “weight sharing” or “symmetry exploitation”—which works here specifically in the linear/quadratic case (forming a “DHPS” (DHP Symmetrized)). But DHPS is identical to GDHP in the linear/quadratic case. GDHP, however, is a general nonlinear adaptation method.

[0301] For DHP, where Ω may be a general positive symmetric matrix that is not necessarily the identity matrix, equation 89 is replaced with: $\begin{array}{cc}\left({\underset{\_}{e}}^T\Omega \frac{\partial }{\partial W_{\mathrm{ij}}}\underset{\_}{e}\right)=\sum _k{\left(\Omega \underset{\_}{e}\right)}_k\frac{\partial }{\partial C_{\mathrm{ij}}}{e}_k& \left(103\right)\end{array}$

[0302] By replacing the vector e on the left side of equation 90 with Ωe, equation 91 changes to:

[0303] <∂_tC>=Ω(V+A^TCA−C)^T  (104)

[0304] The resulting error dynamics are of course:

[0305] <∂_tM>=Ω(A^TMA−M)^T  (105)

[0306] For GDHP, the error dynamics are just the Symmetrized version of equation 105:

[0307] <∂_tM>=Ω(A^TMA−M)^T+^T(A^TMA−M)Ω  (106)

[0308] This provides the framework for the construction of a universal adaptive controller for the MIMO linear deterministic case, guaranteed to stabilize (asymptotically) any unknown plant which is stabilizabie. In other words, this universal adaptive controller should be able to stabilize any plant for which a fixed controller could stabilize (asymptotically) the plant if the dynamics of the plant were known.

[0309] Such an achievement would require total system stability theorems in the spirit of Narendra and Annaswamy [1] and others. To do this requires the individual components—Action, Model and Critic—which are at least stable under the most benign conditions, the conditions where the other components are held fixed. Stable system identifiers and Action/controller systems already exist [1]. However, no one has ever provided Critic adaptation methods which meet the strict standards for stability discussed below.

[0310] In actuality, the present invention focuses on a universal adaptive controller for the linear deterministic case. It utilizes a controller that is based on a Critic adaptation method which converges to the right answer in the stochastic case, derived from an adaptation rule which makes sense in the general nonlinear stochastic situation. Adhering to this constraint should make it easier to carry through the extensions to the stochastic and nonlinear cases in the more distant future.

[0311] Narendra has asserted in several contexts that multiple-model methods are a key to intelligent control. There are several lines of argument which support that viewpoint [13]. Nevertheless, apparently the complexity and the power of multiple models should not be wasted on patching up more elementary problems, like the minimum phase and sign-information problems in adaptive control. When a plant is known to be linear, and when the state space is not divided up into different regions governed by different laws, then single-model methods should suffice.

[0312] There is good reason to believe that the present restrictions on adaptive control (minimum phase and knowledge of gain structure) can be overcome in the linear case, without resorting to multiple models. These problems—as well as slow transient response—are essentially the result of minimizing a simple measure of tracking error at time t+1. An approach which converges to the dynamic programming solution should converge to a stable and optimal solution which, like model-predictive-control should also show good response to transient disturbances, when sufficient input data is available.

[0313] On the other hand, there are three possible uses of multiple models which remain valid and important:

[0314] (1) The use of multiple models to explicitly represent qualitative uncertainty about the world, when that uncertainty cannot be adequately represented by probability distributions over continuous [8,ch.13] or discrete [75] variables. This kind of intellectual strategy is very important to sophisticated symbolic reasoning in humans.

[0315] (2) The use of multiple models to characterize very different regions of state space [13]. Piecewise linear models of nonlinear systems are one classic example of this approach. Piecewise linear models or approximations can emerge in many different ways. For example, Sugeno's classic fuzzy controller for a helicopter (presented at the World Automation Congress 1998) is often cited as the most impressive success story for fuzzy control. It is based on fuzzy IF-THEN rules which translate down, in practice, to a piecewise linear system. Classical gain-scheduling embodies the same underlying principles, with less user-friendly interface. Multiple nonlinear models have also proven very useful and powerful, and are now very popular in parts of the neural network community [76].

[0316] (3) The use of multiple models to characterize different tasks or action schemata or decision blocks [9,13,46] in a hierarchical decision-making system.

[0317] Of these three directions, the application to piecewise linear control seems most relevant to the present state of the art in control engineering. Improved control within each linear region would clearly be very important to these kinds of schemes. By better developing control possibilities within each linear region, it should be easier to build more powerful control systems by linking together these multiple regions.

[0318] The existing literature on stable control for piecewise linear systems (as in [77] or in Sugeno's stability proofs) emphasizes the development of a common quadratic Liapunov function, applicable to all of the regions. When this works, it certainly yields stable controllers. However, it is evident that there are many piecewise linear systems which can be stabilized in principle, in which the same Liapunov function cannot be used across all regions. (For example, one region may be nonminimum phase, relative to the metric or coordinate system implied by the Liapunov function used in another region.) But how piecewise quadratic Liapunov functions be constructed (or even found)? (Or of whatever class of Liapunov function is guaranteed to work in piecewise linear systems that can be stabilized.) The first step in developing such constructive methods is to fully develop constructive methods for individual linear/quadratic regions.

[0319] Critics must interface properly at the boundaries between regions. The new Werbos-Bellman equations given in [13] (extracted from a longer, denser paper [46]) provide a rigorous basis for managing that interface. The development of those equations was initially motivated by an effort to develop more truly brain-like intelligent control: however, the mathematics is not limited to that application.

[0320] At least six areas touch on the present invention:

[0321] (1) Even in the “simple” testbed problem of supervised learning (learning a static mapping from an input vector X to a vector of targets or dependent variables Y), improved learning speed and generalization ability, exploiting concepts such as “syncretism” and “simultaneous recurrence” [32] are needed. Also see [34].

[0322] (2) In adapting Model networks, architectures (like SEDP [8,ch. 13] and like [75]) which can represent probability distributions, not just deterministic input-output relations, need to be explored. “Robustness over time” [8,ch. 10] needs to be improved as well. Error Critics [8, ch. 13:26:15] need to be explored as an alternative to BTT, in order to combine real-time learning with the robust kind of memory now available only with BTT. The improvements needed in supervised learning should be extended to Model training as well.

[0323] (3) Substantial improvements in basic ADP technology are also needed.

[0324] (4) At a higher level, designs for systems that perform “temporal chunking” in an adaptive way need to be developed, understood, tested and applied [9,13,46].

[0325] (5) Better systems are needed for a “brain-like” stochastic search, in order to help decision-making system escape from local minima, without the artificial limitations of methods like genetic algorithms that do not appear to be truly brain-like. For example, when trying to find actions or decisions u which maximize U(α,u) for a class of static function maximization problems, in which each specific problem is characterized by set of parameters α, how should the problem be attacked? (As an example, α might represent the locations of cities in a traveling salesman problem.) An Option Network [46] can be trained so that it inputs α and outputs possible u, in a stochastic way, and tries to ensure that it outputs u according to a Gibbs distribution, c*exp(−(k/T)U), where T is a temperature parameter also input to the network.

[0326] (6) Principles like spatial symmetry and objective symmetry [46] need to be implemented in new classes of neural networks, in order to make it possible for these networks to cope with the huge volumes of input which biological organisms can handle. In effect, neural adaptive designs need to learn to do what some A1 systems already do in a crude, brute-force manner; build up spatial “world models;” analyze objects one at a time, in a multiplexing kind of mode; perform spatial chunking; and input and output networks of relations rather than fixed-length vectors. Again, see [46].

[0327] The effort to really understand the level of intelligence found in basic mammalian brains, in a mathematically sophisticated way, is already a difficult enough challenge for now. It may be a prerequisite to a deeper, less fuzzy understanding of the higher levels.

[0328] To produce the most robust system, it would have been preferable that both of the following basic test were met:

[0329] (1) “Correct equilibrium”, defined as follows: when C equals the correct C* matrix, defined by equation 56, then the adaptation methods results in <∂_tC>=0 in general in the linear stochastic case.

[0330] (2) Stable convergence to the equilibrium: The adaptation rule should possess a property called quadratic unconditional stability (to be discussed below) in the linear deterministic case, at least.

[0331] Unfortunately, none of these methods passed both of these tests. In particular, none of the Galerkinized methods discussed above possess the property of correct equilibrium. (The results of HDPC and HDPGC are not known precisely.) These results are consistent with earlier work [17], addressing the case of linear/lineal multivariate plants; in that work HDP demonstrated correct equilibrium for stochastic plants, but HDPG did not.

[0332] From the definition of M in equation 63, the property of “correct equilibrium” as defined above is equivalent to the property that <∂_tM>=0 when M=0. For HDP, the dynamics of M in the stochastic case are given in equation 83. When M=0 is substituted into equation 83, the right-hand side equals zero. Thus HDP possesses correct equilibrium. Likewise, the dynamics of M in the stochastic case for DHP are given in equation 97. Again, inserting M=0 sends the right-hand side to zero, which demonstrates correct equilibrium.

[0333] Finally, for GDHP, in the most general version analyzed above, the dynamics of M in the stationary case are given in equation 105. Once again, when M=0<∂_tM>=0, as required. Thus GDHP also possesses correct equilibrium.

[0334] The equilibrium properties of these methods can again be analyzed for HDPG, DHPG, and GDHPG in the same way as for GDHP. The dynamics of M for HDPG in the stochastic case are given in equation 84. Substituting M=0 into that equation produces:

[0335] <∂_tM>=−2QC*Axx^T−2Axx^TA^TC*Q+Tr(C*Q)Q−<(w^TC*w)ww^T>  (107)

[0336] If Q=0 (the deterministic case), this will still go to zero. This is like the results in [17], where it was shown that HDPG yields an incorrect equilibrium almost always. There was a small list of special cases, of plants where the equilibrium would be correct. At the top of that list was the case where Q=0. But beyond that, for normal plants C* is not expected to be a singular matrix. In a stochastic system, it is not expected that the state vectors x will always fit the condition Ax=0 or (even less likely) always have the property that they exactly cancel out the last two terms in equation 107. It is clear that equation 107 will usually (in fact, almost always) be different from zero for a normal stochastic system. Certainly equation 107 is not always zero, as the property of “correct equilibrium” requires.

[0337] Likewise, the dynamics of M for DHPG in the stochastic case are given in equation 98. Substituting M=0 into that equation produces:

[0338] <∂_tM>=−AA^TC*Q  (108)

[0339] Once again, when Q is nonzero, this matrix is expected to be nonzero as well. Note that this equation indicates a constant, unremitting bias away from the correct equilibrium, regardless of the states x actually experienced.

[0340] For GDHPG (with ω=I), the error dynamics for the stochastic case are given in equation 102. Substituting M=0 yields:

[0341] <∂_tM>=−AA^TC*Q−QC*AA^T  (109)

[0342] The conclusions here are exactly the same as for DHP. In summary, none of the Galerkinized methods discussed above possess “correct equilibrium” in the linear stochastic case.

[0343] Which methods then possess quadratic unconditional stability in the deterministic case? Preliminarily, the concept of stability as used herein is extremely stringent—far more stringent than needed in practice. However, (1) it is a relatively easy point of departure for an initial phase of research; and (2) methods that do meet this stringent test tend to be of practical interest as well.

[0344] There are many, many concepts of stability useful in different types of application [1,73, 78]. To define the type of quality discussed herein, the term “q-stability” is used. Systems that do not possess q-stability will be called “q-unstable.”

[0345] Q-stability is defined as follows, for the kinds of adaptive systems discussed above. For any dynamic process {A,Q,V}, an adaptation process for C is q-stable if there exists a Liapunov function)—or, equivalently, Λ(M,A,Q,V)—such that:

[0346] (1) Λ is a positive definite quadratic function of M;

[0347] (2) <∂_tΛ>≦0 for all combinations of M and x.  (110)

[0348] In equation 110, the angle brackets refer to the expectation value with respect to the random vector w.

[0349] Intuitively, unconditional stability requires stability in the face of sequence of inputs x—not just the sequences that the plant would normally generate if it obeys equations 45 or 46 (or even 47 or 48). In effect, the system is to learn and thrive even when the experience it grows up on is “a movie made by a devil” —a sequence of collages (pairs of x(t) and x(t+1)), each of which flows naturally internally, but which fade in and out in an illogical manner, specifically chosen so as to destabilize the learning process. Q-stability is the special case of unconditional stability, in which the Liapunov function is quadratic. (An example is given later of a linear dynamic system which is unconditionally stable but not q-stable.)

[0350] The notion of q-stability given above is almost identical to the notion of quadratic stability as discussed in the seminal book by Boyd et al [79]. On page 73 of that book, they note that it is possible for a Linear Differential Inclusion (like the dynamic rules for M) to be stable, without being quadratically stable. On page 74, they cite a number of examples; they state that unconditional stability in a more general sense occurs “. . . if and only if there is a convex Liapunov function that proves it . . . IN general, computing such Liapunov functions is computationally intensive, if not intractable.”

[0351] Ideally, adaptive methods should be found that are q-stable for all dynamic processes such that the eigenvalues of A lie within the unit circle. The Galerkinized methods do possess that property, for deterministic linear plants. But the other methods—the methods which possess “correct equilibrium” —are q-unstable, even in the deterministic case, except for a “set of measure zero” in the universe of dynamical systems. For practical purposes, they are q-stable only when A is a perfectly symmetric matrix.

[0352] Q-stability applies to a well-known, fundamental learning rule used in supervised learning and signal processing. That rule is the real-time learning rule called Least Mean Squares (LMS) developed by Widrow and Hoff [80], in a linear multivariate case, which is known to have robust stability.

[0353] Assume that a sequence of pairs, x(t) and y(t), of a vector x and a scalar y are observed that are governed by the stochastic process:

[0354] y(t)=(b*)^Tx(t)+w(t),  (111)

[0355] where w is a normally distributed random variable. Also assume a learning process based on: $\begin{array}{cc}\Delta b_i=-\alpha \frac{\partial E}{\partial b_i},& \left(112\right)\end{array}$

[0356] where:

[0357] E=½e²=½(y−b^Tx)²  (113)

[0358] Substituting 113 into 112, produces: $\begin{array}{cc}\frac{\partial E}{\partial b_i}=-\left(y-{\underset{\_}{b}}^T\underset{\_}{x}\right)x_i& \left(114\right)\end{array}$

[0359] Assuming infinitely small learning rates, equation 114 is substituted into equation 112, and represents the result as a differential equation:

[0360] ∂_tb=(y−b^Tx)x  (115)

[0361] Substituting equation 111 into equation 115 produces:

[0362] ∂_tb=(b^Tx+w−b^Tx)x  (116)

[0363] In parallel with equation 63, a new vector δ: may be implicitly defined that represents the error of the adaptive system in estimating the parameter vector b:

[0364] b=b*+δ  (117)

[0365] Substituting this definition into equation 116, and taking expectation values with respect to the random variable w, produces:

[0366] <∂_tδ>=−(δ^Tx)x  (118)

[0367] The following Liapunov function is then selected for this system:

[0368] Λ(δ)=½|δ|²=½δ^Tδ  (119)

[0369] This meets the first condition for a Liapunov function, that it is positive δ definite in 8 (which plays the same role as M above). To verify the second condition, equations 118 and 119 are used to compute:

[0370] <∂_tΛ>=<δ^T∂_tδ>=δ^T<∂_tδ>=−(δ^Tx)(δ^Tx)≦0  (120)

[0371] For all four systems, stability occurs using the same Liapunov function:

[0372] Λ(M)=½Tr(M²)  (121)

[0373] Because M is a symmetric real matrix, its square is guaranteed to be a symmetric non-negative matrix. All of its eigenvalues are guaranteed to be nonnegative. Λ is positive definite, because the trace—the sum of the eigenvalues—can only be zero if all of the eigenvalues are zero, which would imply M=0. Therefore, equation 121 meets the first condition for a Liapunov function. To complete the proof for each of the four Galerkinized methods to be discussed here, only the second condition need be verified for each of them.

[0374] Because M is a symmetric matrix, differentiating equation 121 produces:

[0375] ∂_tΛ(M²)=Tr(M(∂_tM))  (122)

[0376] The dynamics of M, the error matrix for the weights, in HDPGC are given in equation 66. Substituting that equation into equation 122 (and choosing α=1 for simplicity), produces:

[0377] δ_tΛ=−Tr(Me(xx^TA^T+Axx^T))=−e(Tr(Mxx^TA^T)+Tr(MAxx^T))  (123)

[0378] Exploiting the well known theorem (easily proved) that Tr(AB)=Tr(BA), even for nonsquare matrices A and B, (123) reduces to:

[0379] ∂_tΛ=−e(x^TA^TMx+x^TMAx)  (124)

[0380] Substituting in from equation 64, this is simply:

[0381] ∂_tΛ=−(x^TA^TMx+x^TMAx)²≦0  (125)

[0382] This completes the proof of u-stability.

[0383] The error dynamics for HDPG in the stochastic case are given in equation 84. For the deterministic case, where Q=0 and w=, equation 84 reduces to:

[0384] ∂_tM=d(xx^T−Axx^TA^T),  (126)

[0385] where the following new scalar is defined in order to simplify the algebra:

[0386] d=x^T(A^TMA−M)x  (127)

[0387] Substituting equation 126 into equation 122 produces:

[0388] ∂_tΛ=Tr(M(∂_tM))=d Tr(Mxx^T−MAxx^TA^T)=dx^TMx−dx^TA^TMAx=dx^T(M−A^TMA)x=−d²≦0  (128)

[0389] Again, QED.

[0390] The error dynamics for DHPG in the stochastic case are given in equation 98. For the deterministic case, where Q=0, equation 98 reduces to:

[0391] ∂_tM=(A^TMA−M)xx^T−A(A^TMA−M)xx^TA^T  (129)

[0392] If, for convenience, the symmetric matrix is defined as:

[0393] D=A^TMA−M,  (130)

[0394] equation 129 is simplified to:

[0395] ∂_tM=Dxx^T−ADxx^TA^T  (131)

[0396] Following equation 122, the following relationship derived:

[0397] ∂_tΛ=Tr(M(∂_tM))=Tr(MDxx^T)−Tr(MADxx^TA^T)=x^TMDx−x^TA^TMADx=x^T(M−A^TMA)Dx=−x^TD²x=−|Dx|²≦0  (132)

[0398] QED.

[0399] The error dynamics for GDHPG in the stochastic case, for ω₀=0 and ω=1, are given in equation 102. Equation 102 is really just the symmetrized version of equation 98. In the deterministic case, where Q=0, it reduces to the symmetrized version of equation 129, which can be expressed as follows:

[0400] M=Dxx^T−ADxx^TA^T+xx^TD−Axx^TDA^T,  (133)

[0401] using the definition of “D” given in equation 130. Following equation 122, the following relationship is derived:

[0402] ∂_tΛ=Tr(M(∂_tM))=Tr(MDxx^T)−Tr(MADxx^TA^T)+Tr(Mxx^TD)−Tr(MAxx^TDA^T)=x^TMDx−x^TA^TMADx+x^TDMx−x^TDA^TMAx=x^T(M−A^TMA)Dx+x^TD(M−A^TMA)x=−x^TD²x−x^TD²x=−2|Dx|²≦0  (134)

[0403] Once a Liapunov function, has been found, q-stability can be proven. For a given plant, instability can also be proven in a straightforward way, by giving examples of training schedules {x(t)} which cause an explosion in the estimation errors M. However, to prove q-instability across a broad range of possible plants is more difficult. It requires a more abstract mathematical approach.

[0404] First of all, our focus here will be on q-instability for four methods—HDPC, HDP, DHP and GDHP—in the deterministic linear case. This reduces down to an analysis of four different dynamical systems for the matrix M—equation 65 (for HDPC), equation 73 (for HDP), equation 105 (with angle brackets removed, for generalized DHP) and equation 106 likewise, for generalized GDHP). It can easily be seen that these dynamical systems are linear in the matrix M. Thus the easiest way to analyze these systems is by applying linear systems theory to them. In this case, it requires a certain amount of care to do this.

[0405] The original dynamical system (equations 44 or 45) involved a state vector, x, with n components. In formal terms, xεRⁿ. But the linear deterministic systems of equations 65,7 73, 105 and 106 concern a matrix M. In order to apply ordinary linear systems theory, M, in equation 105, is treated as a kind of vector of dimensionality N=n². In other words, MεR^Nis the representation of the matrix M as a vector in n²-dimensional space. For equations 65, 73 and 106, M is taken from the space of symmetric matrices, which only have N′=n(n+1)/2 independent degrees of freedom. Then in those cases, the vector representation of M will be M εR^N′. Throughout, the quantities n, N′ and N are defined in this way. It is also noted that in this formulation, the usual matrix Ricatti equations are ordinary linear equations for the vector C in R^N′. This permits exact finite-time solutions using ordinary methods for the solution of linear equations over R^N′.

[0406] For q-stability in the linear case, Λ(C,Q,V,A) is required to be quadratic in C. The Liapunov function can be expressed as a quadratic about the desired stable equilibrium, C*; in other words, it can be expressed as a quadratic in M, if such a Liapunov function exists. Because Q=0 in the linear deterministic case, and because the dynamics of M in equations 65, 73. 105 and 106 do not depend upon V, this means that Λ(M,A) is required to be quadratic in M. Thus a Liapunov function in the following form is used: $\begin{array}{cc}A\left(M,A\right)=\frac{1}{2}\sum _{\mathrm{ijkl}}{T}^{\mathrm{ijkl}}\left(A\right)M_{\mathrm{ij}}{M}_{\mathrm{kj}}& \left(135\right)\end{array}$

[0407] (Even though T^ijklshould always be understood to be a function of A, sometimes this dependence is only implicit in the description below.)

[0408] Without loss of generality, T can be required to be symmetric in M. In other words:

[0409] T^ijkl=T^kilj, for all i,j,k,l  (136)

[0410] Furthermore, if we think of M as a vector in R^Nor R^N′, we may also define:

[0411] $\begin{array}{cc}\underset{\_}{L}=\nabla \Lambda =\sum _{i,j}{T}^{\mathrm{ijkl}}\left(A\right)M_{\mathrm{ij}},& \left(137\right)\end{array}$

[0412] where L also corresponds to a matrix L. In the case where M is a vector in R^N′, without a loss of generality, it may be required that:

[0413] T^ijkl=T^jikl=T^ijlk=T^jilk,  (138)

[0414] and that L is therefore a symmetric matrix.

[0415] The first condition for q-stability requires that T(A) must be positive definite, when considered as an N-by-N or N′-by-N′ matrix. The second condition requires that: $\begin{array}{cc}{\partial }_t\Lambda =\sum _{\mathrm{ijkl}}{T}^{\mathrm{ijkl}}{M}_{\mathrm{ij}}\left({\partial }_t{M}_{\mathrm{kl}}\right)=\mathrm{Tr}\left(L^T{\partial }_{t}M\right)\le 0& \left(139\right)\end{array}$

[0416] for all M and x, for the particular matrix A under study. Note that equation 139 specifically follows for the choice of T as a symmetric matrix. It is not possible to satisfy the requirement for q-stability with a matrix T which is not symmetric, unless it is also possible to satisfy it for the symmetrized version of that matrix. Thus to prove q-instability in general, it is enough to prove that no such symmetric T exists.

[0417] Could the requirement in equation 138 be satisfied for some symmetric T, in HDPC? To evaluate this, substitute equation 65 into equation 138 (with α=½), which yields the requirement:

[0418] Tr(L(½exx^T))≦0 for all x, M  (140)

[0419] Note that L is a symmetric matrix, in this case, because M is a symmetric matrix, and because equation 138 is imposed.

[0420] Substituting equation 64 into 140, the requirement of equation 138 reduces to:

[0421] (x^TLx)(x^T(A^TM+MA)x)≦0 for all x and M  (141)

[0422] The analysis of equation 141 leads to the following lemma:

[0423] Lemma on L and B: If:

[0424] (x^TLx)(x^TBx)≦0 for all x  (142)

[0425] for some symmetric matrices L and B, then either: (1) B is positive semidefinite and L is negative semidefinite: (2) L is positive semidefinite and B is positive semidefinite; or (3) L=−kB for some nonnegative scalar k.

[0426] Proof: Since L is real and symmetric, without a loss of generality, it may be assumed that a similarity transformation is performed on L, B and x based on the eigenvectors of L, such that L becomes a diagonal matrix. After that, if the elements on the diagonal of L are all nonnegative or all nonpositive, then x^TLx is always negative or always positive. In those cases, equation 142 requires that x^TBx is always nonpositive or always nonnegative, respectively, for all x, which then gives us cases (2) and (1) respectively. The only alternative possibility is that L contains a mix of diagonal elements (eigenvalues), some positive and some negative.

[0427] Let e_iand e_jbe the basis vectors (after the similarity transformation) corresponding to one of the positive eigenvalues λ_iof L, and one of the negative eigenvalues, λ_j, respectively. When x=e_i, equation 142 states that B_ji≦0. Likewise, when x=e_jthen B_ji≧0. For x=e_i+ye_j, it can be shown that:

[0428] x^TLx=λ_i+y²λ_j  (143a)

[0429] x^TBx=B_ii+2yB_ij+y²B_jj  (143b)

[0430] The graph for equation 143a as a function of y is a parabola, crossing the x-axis at two points, where y=±sqrt(−λ_i/λ_j). Equation 142 shows that the quadratic expression in 143b must cross the x-axis at the same two points. From elementary algebra, this means that B_ij=0 and B_ii/B_jj=λ_i/λ_j. Applying this argument to all possible positive/negative λ_i/λ_jpairs (including λ_k/λ_jinvolving different k but the same j), it may be deduced that the entire diagonal of B must equal −k multiplied by L, for some nonnegative scalar k. It may also be deduced that all the B_ijcrossterms for λ_ipositive and λ_jnegative must be zero.

[0431] Next, consider x of the form e_i+ye_j+ze_k, where λ_k, λ_i>0 and λ_j<0. This produces two expressions quadratic in y, analogous to 143a and 143b, which must match crossing points as before, for every possible value of z. The only way this can happen is if B_ik=0. In a similar way, it is known that the offdiagonal terms in the negative block of B are also zero. In summary, all of the offdiagonal terms of B must be zero—the ones which cross between the negative and positive blocks, and the ones which go between different elements within each block. Thus B must be a diagonal matrix. Since the diagonal of B must equal −kL, this completes the proof of the lemma.

[0432] Applying the above lemma to equation 141, it is required that:

[0433] L(A,M)=−k(A,M)(A^TM+MA),  (144)

[0434] for all matrices M such that A^TM+MA is not positive semidefinite or negative semidefinite.

[0435] Assume MεP(A) when A^TM+MA is positive semidefinite. It is straightforward to prove that P(A) is a convex cone[81] by simply applying the definition of a convex cone: (1) if MεP(A), then kMεP(A) (for k≧0) because k multiplied by a positive semidefinite matrix is still a positive semidefinite matrix: (2) if M₁εP(A) and M₂εP(A), then (I−k)M₁+kM₂(for 1>k>0) is in P(A), because the sum of two positive definite matrices is also positive definite. Also, MεN(A) iff −MεP(A); this is also a closed convex cone. If K(A) is defined as the rest of R^a, after P(A) and N(A) are removed; clearly it, too, is a connected cone.

[0436] Equation 137 shows that L would be required to be linear in M. Consider any two matrices, M₁and M₂, both in K(A), such that M₁+M₂is also in K(A), and that A^TM₁+M₁A is not just equal to a scalar multiple of A^TM₂+M₂A. In this case, the linearity of L as a function of M requires that k(A,M₁)=k(A,M₂)=k(A,M₁+M₂). Likewise, the linearity of L requires that k(A,cM)=k(A,M) for any scalar c. Since all matrices M in K(A) are connected to each other by some combination of adding and scaling. k(A,M) is actually k(A)—a constant with respect to M, all across K(A). Furthermore, all matrices in P(A) or N(A) can be reached as linear combinations of matrices in K(A). (In other words, it is easy to write any positive definite symmetric matrix as a linear combination of two symmetric matrices which are not definite. A convex combination is not required, only a linear combination.) Therefore:

[0437] L(A,M)=−k(A)(A^TM+MA) for all M,  (145)

[0438] for any linear dynamic process (A) which is q-stable.

[0439] Combining equations 137 and 145, produces the requirement that: $\begin{array}{cc}\sum _{\mathrm{kl}}{T}^{\mathrm{ijkl}}\left(A\right)M_{\mathrm{kl}}=-k\left(A\right)\sum _k\left(A_{\mathrm{ik}}^T{M}_{\mathrm{kj}}+M_{\mathrm{ik}}{A}_{\mathrm{kj}}\right)\mathrm{for}\mathrm{all}M& \left(146\right)\end{array}$

[0440] Differentiating with respect to M_kl(recalling that M_kl=M_lkis one element of MεR^N′) and exploiting the fact that T can be multiplied by a positive scalar without changing its validity as a Liapunov function. produces the requirement that:

[0441] T^ijkl(A)=−(δ_jlA_ik^T+δ_jkA_il^T+δ_ilA_kj+δ_ikA_lj)  (147)

[0442] Again, equation 147 has been deduced as a necessary requirement for a symmetric, quadratic Liapunov function T for a deterministic linear system. But equation 136 is another necessary requirement. The combination of these two requirements, which apply to all combinations of i,j,k and l, implies a requirement that:

[0443] δ_jlA_ik^T+δ_jkA_il^T+δ_ilA_kj+δ_ikA_lj=δ_ljA_ki^T+δ_liA_kj^T+δ_kjA_il+δ_kiA_jl  (148)

[0444] For example, in the case where i=j=k, this requires that:

[0445] δ_ilA_ii^T+A_il^T+δ_ilA_ii+A_li=δ_liA_ii^T+δ_liA_ii^T+A_il+A_il  (149)

[0446] which reduces to:

[0447] 2A_li=2A_ilfor all i and l  (150)

[0448] which is satisfied if and only if the matrix A is symmetric.

[0449] The conclusion of this analysis is that HDPC can be q-stable only in the case where A happens to be symmetric. Matrices of this sort are a set of measure zero in the space of all matrices A! More importantly, there are few real-world applications where the dynamic matrix A is exactly symmetric. This is far more restrictive than the restrictions described for the old IAC designs. When A does happen to be symmetric, it can be shown that equation 147 gives a valid Liapunov function. HDPC clearly fails to meet other goals described herein, so further discussion of it is omitted.

[0450] Notice that L, the gradient of Λ, is defined for any differentiable Liapunov function over M. Thus equations 140, 141 and 144 are valid requirements for unconditional stability in general, not just q-stability. Q-stability was used in this argument only to derive the narrower requirement given in equation 145. Prior to the development of this proof, a variety of two-by-two and diagonalized (complex) examples were analyzed, where the results were all consistent with the general result given here. Intuitively, the argument here is saying that the required gradient field L is not an integrable field, because it does not meet the usual condition for the integrability of a vector field. In other words, the required L cannot actually be the gradient of a scalar field Λ, as required, because its derivatives do not possess the required symmetry.

[0451] Equation 144 seems very close to equation 145, on the surface; however, this may be misleading. There appear to be some dynamical systems A, which are not exactly symmetric, where HDPC is unconditionally stable, based on a Liapunov function which is not quadratic. However, this does not imply unconditional stability across all controllable systems A.

[0452] The proof of q-instability of HDP proper follows exactly the same structure as the earlier proof. After all, the dynamics of M in equation 73 are virtually identical to the dynamics in equation 65, except the definition of “e” changes (from equation 72 instead of 64).

[0453] More precisely, the proof of q-instability for HDP is exactly the same except that “A^TM+MA” is replaced everywhere by “A^TMA−M”. In place of equation 146, it is required that: $\begin{array}{cc}\sum _{\mathrm{kl}}{T}^{\mathrm{ijkl}}{M}_{\mathrm{kl}}=-k\left(A\right)\left(-M_{\mathrm{ij}}+\sum _{\mathrm{kl}}{A}_{\mathrm{ik}}^T{M}_{\mathrm{kl}}{A}_{\mathrm{lj}}\right)& \left(151\right)\end{array}$

[0454] for all M. Once again, differentiating with respect to M_kl=M_lkand exploiting the invariance of T with respect to positive scalar multiplication, derives the requirement:

[0455] T^ijkl=δ_ikδ_jl+δ_ilδ_jk−(A_ik^TA_lj+A_il^TA_kj)  (152)

[0456] Once again, equation 134 then requires:

[0457] δ_ikδ_jl+δ_ilδ_jk−A_ik^TA_lj−A_il^TA_kj=δ_kiδ_lj+δ_kjδ_li−A_ki^TA_jl−A_kj^TA_il  (153)

[0458] for all combinations of i, j, k, and l. For example, choosing k=l and summing the results over k, equation 15e reduces to the necessary requirement that:

[0459] −A^TA−A^TA=−AA^T−AA^T  (154)

[0460] Thus A commutes with A^T, and both of them commute with A+A^T. Thus, A is a “normal matrix”[81,82] its eigenvectors (like those of A+A^T) will usually be real and orthogonal. (Exceptions may be possible if A+A_Thas degenerate eigenvalues.) When a real matrix has eigenvectors which are all real and orthogonal, it must be symmetric. Once again, the matrix A must either be symmetric, or it must fall into some other set of exceptions which is also of measure zero.

[0461] The requirement given in equation 154 is only a necessary requirement for q-stability, not a sufficient requirement. The necessary requirement is already strong enough to show that HDP does not provide q-stability for a broad enough range of dynamical systems.

[0462] Because GDHP preserves the symmetry of M, while DHP does not, it is easier to begin by considering GDHP as a next system to be examined. For GDHP, q-stability again requires a quadratic Liapunov function that satisfies equations 136 through 139. In the deterministic case, the error dynamics for GDHP as given in equation 106 may be written as:

[0463] ∂_tM=Ω(A^TMA−M)xx^T+xx^T(A^TMA−M)Ω=Dxx^T+xx^TD^T.  (155)

[0464] where define:

[0465] D=Ω(A^TMA−M)  (156)

[0466] Substituting equation 155 into equation 139, and recalling that L must be symmetric in this case, gives the requirement that:

[0467] Tr(LDxx^T+Lxx^TD^T)=x^TLDx+x^TD^TLx=2(Lx)^T(Dx)≦0  (157)

[0468] for all vectors x.

[0469] It is safe at this point to assume that D and L are both nonsingular, for the following reasons. First, the case A=±I is less important since the goal is to show that GDHP is almost always q-unstable, when A is not symmetric. Second, for a sensible adaptation algorithm, Ω would be required to be nonsingular. Third, Λ is required to be positive definite, which implies that L must be nonsingular when M is. In summary, for all nonsingular M, L and D must be nonsingular (for q-stability), produces the matrix:

[0470] D=S(A,M)L  (158)

[0471] Inserting this definition into equation 157 shows that S(A,M) must be a negative semidefinite matrix, for all M and A. (By continuity, this is expected for all matrices M.)

[0472] Combining equation 158 with equations 137 and 156, the requirement for q-stability becomes: $\begin{array}{cc}{L}_{\mathrm{ij}}=\sum _{\mathrm{kl}}{T}^{\mathrm{ijkl}}{M}_{\mathrm{kl}}={\left(S^{-1}\Omega \left(A^T\mathrm{MA}-M\right)\right)}_{\mathrm{ij}}& \left(159\right)\end{array}$

[0473] The analysis of equation 159 leads to the following lemma of tensor linearity. Suppose that: $\begin{array}{cc}\sum _{\mathrm{kl}}{P}^{\mathrm{ijkl}}{M}_{\mathrm{kl}}=\sum _{\mathrm{klv}}{B}^{\mathrm{ivkl}}{M}_{\mathrm{kl}}{C}_{\mathrm{vj}}\left(M\right)& \left(160\right)\end{array}$

[0474] for all symmetric matrices M, for some tensors P and B which are not functions of M. Assume also that B has the property that BM is nonsingular almost everywhere, in the space of MεR^N′. Then C_vj(M)=C_vj(0) is not actually a function of M.

[0475] Proof: The validity of equation 160 implies the validity of equation 160 differentiated with respect to M_αβ, for any α and β; this yields: $\begin{array}{cc}{A}^{\mathrm{ij}\mathrm{\alpha \beta }}=\sum _v{B}^{\mathrm{iv}\mathrm{\alpha \beta }}{C}_{\mathrm{vj}}\left(M\right)+\sum _{\mathrm{klv}}{B}^{\mathrm{ivkl}}{M}_{\mathrm{kl}}\frac{\partial C_{\mathrm{vj}}}{\partial M_{\mathrm{\alpha \beta }}}& \left(161\right)\end{array}$

[0476] Multiplying by M_αβ and summing over α and β, produces the requirement that: $\begin{array}{cc}\sum _{\mathrm{\alpha \beta }}{A}^{\mathrm{ij}\mathrm{\alpha \beta }}{M}_{\mathrm{\alpha \beta }}=\sum _{\mathrm{\alpha \beta }v}{B}^{\mathrm{iv}\mathrm{\alpha \beta }}{M}_{\mathrm{\alpha \beta }}{C}_{\mathrm{vj}}\left(M\right)+\sum _{\mathrm{klv}\mathrm{\alpha \beta }}{B}^{\mathrm{ivkl}}{M}_{\mathrm{\alpha \beta }}{M}_{\mathrm{kl}}\frac{\partial C_{\mathrm{vj}}}{\partial M_{\mathrm{\alpha \beta }}}& \left(162\right)\end{array}$

[0477] But the first two terms in equation 162 ire essentially identical to equation 160. Subtracting out equation 160 from equation 162, produces: $\begin{array}{cc}\sum _{\mathrm{klv}\mathrm{\alpha \beta }}{B}^{\mathrm{ivkl}}{M}_{\mathrm{\alpha \beta }}{M}_{\mathrm{kl}}\frac{\partial C_{\mathrm{vj}}}{\partial M_{\mathrm{\alpha \beta }}}=0& \left(163\right)\end{array}$

[0478] Defining two new matrices produces: $\begin{array}{cc}{D}_{\mathrm{iv}}=\sum _{\mathrm{kl}}{B}^{\mathrm{ivkl}}{M}_{\mathrm{kl}}& \text{(164a)}\\ R_{\mathrm{vj}}=\sum _{\mathrm{\alpha \beta }}{M}^{\mathrm{\alpha \beta }}\frac{\partial C_{\mathrm{vj}}}{\partial M_{\mathrm{\alpha \beta }}}& \text{(164b)}\end{array}$

[0479] Using these definitions, equation 103 may be expressed as:

[0480] D(M)R(M)=0  (165)

[0481] Because BM is assumed to be nonsingular almost everywhere in M, D(M) will be nonsingular almost everywhere. Thus almost everywhere in M, equation 165 requires:

[0482] R(M)=0  (166)

[0483] By the requirement for continuity with respect to M, R(M)=0 is required for all M. Equatione 166, plus the definition in equation 164b, are an instance, in R^N′, of the equations:

[0484] x^T∇f(x)=0 for all x  (167)

[0485] Integrating equation 167 along any line radiating out from the origin to a point x implies f(x)=f(0). In other words, f is actually a constant with respect to x. Applying this to equations 164b and 166, C_vj(M)=C_vj(0) is not actually a function of M.

[0486] Applying equation 160 and the lemma above to the analysis of equation 159, it may be deduced that S⁻¹cannot actually be a function of M. However, is BM nonsingular almost everywhere in this case? In this case, BM is simply the right-hand side of equation 159, which, for q-stability, must equal L. But in the discussion leading tip to equation 158, it was explained why L is required to be nonsingular almost everywhere in M.

[0487] Accepting that S⁻¹cannot be a function of M, a new matrix is defined as:

[0488] R(A)=S⁻¹(A)Ω  (168)

[0489] Inserting this definition into equation 159, produces:

[0490] L=R(A)(A^TMA−M),  (169)

[0491] but L is required to be symmetric for all symmetric matrices M. Assuming that the N′-by-N′ #matrix which multiplies M to generate A^TMA−M is not singular, this implies that A^TMA−M ranges over all possible symmetric matrices as well (again, excluding consideration of the case A=±1) If R(A)D is symmetric for all symmetric matrices D, then R(A) must be of the form kl, where k is a scalar. Substituting this back into equation 169, and recalling that T is only defined up to a positive scalar anyway (and using some knowledge about systems analysis to set the sign), produces the requirement that: $\begin{array}{cc}\sum _{\mathrm{kl}}{T}^{\mathrm{ijkl}}{M}_{\mathrm{kl}}=M_{\mathrm{ij}}-\sum _{\mathrm{kl}}{A}_{\mathrm{ik}}^T{M}_{\mathrm{kl}}{A}_{\mathrm{lj}}& \left(170\right)\end{array}$

[0492] Differentiating with respect to M_kl=M_lk, produces the requirement:

[0493] T^ijkl=δ_ikδ_jl−A_ik^TA_lj+δ_ilδ_jk−A_il^TA_kj  (171)

[0494] Once again, equation 136 is used to deduce the requirement:

[0495] δ_ikδ_jl−A_ik^TA_lj+δ_ilδ_jk−A_il^TA_kj=δ_kiδ_lj−A_ki^TA_jl+δ_kjδ_li−A_kj^TA_il  (172)

[0496] for all i, j, k and l. This reduces to:

[0497] A_ik^TA_lj+A_il^TA_kj=A_ki^TA_jl+A_kj^TA_il  (173)

[0498] Consider the case k=l and sum over 1 to deduce the necessary requirement:

[0499] A^TA+A^TA=AA^T+AA^T  (174)

[0500] This is essentially the same as equation 154. Q-stability requires that A must either be symmetric, or must belong to some other set of exceptional cases of measure zero. Thus GDHP, like HDPC and HDP proper, is q-unstable for almost all matrices A, under the conditions described herein.

[0501] For q-stability of DHP, equations 135-137 and 139 would still need to hold; however M and L are no longer required to be symmetric, and equation 138 does not apply. Substituting the adaptation rule for generalized DHP in the deterministic case (equation 105 with angle brackets removed) into equation 139, produces the requirement for q-stability that:

[0502] Tr(L^T∂_tM)=Tr(L^TDxx^T)=x^TL^TDx=(Lx)^T(Dx)≦0  (175)

[0503] where D is again defined as in equation 156. Note the similarity of the final inequality in equation 175 to the final inequality in equation 157. Starting from that inequality, the exact same logic is used as before, in order to deduce that equations 158 and 159 are also required here. Likewise, the lemma of tensor linearity can be applied. Thus, q-stability for DHP will require that equation 159 hold, for S⁻¹not a function of M. This again allows the definition in equation 168, which leads to equation 169.

[0504] For DHP, however, it is not possible to deduce that R(A) must equal kl. The requirement for q-stability given in equation 169 may be written out more exhaustively and combined with equation 137, to yield the requirement: $\begin{array}{cc}{L}_{\mathrm{ij}}=\sum _{\mathrm{ijv}}{R}_{\mathrm{kv}}{A}_{\mathrm{vi}}^T{M}_{\mathrm{ij}}{A}_{\mathrm{jl}}-\sum _v{R}_{\mathrm{kv}}{M}_{\mathrm{vl}}=\sum _{\mathrm{ij}}{T}^{\mathrm{ijkl}}{M}_{\mathrm{ij}}\mathrm{for}\mathrm{all}M& \left(176\right)\end{array}$

[0505] Differentiating with respect to M_ij, produces the requirement: $\begin{array}{cc}\sum _v{R}_{\mathrm{kv}}{A}_{\mathrm{vi}}^T{A}_{\mathrm{jl}}-R_{\mathrm{ki}}{\delta }_{\mathrm{lj}}=T^{\mathrm{ijkl}}& \left(177\right)\end{array}$

[0506] Applying equation 136, produces the requirement, for all i, j, k and l, that:

[0507] (RA^T)_kiA_jl−R_kiδ_ij=(RA^T)_ikA_lj−R_ikδ_jl  (178)

[0508] Consider the case where lj. If all of the crossterms A_ljfor lj are zero, the A would be a symmetric matrix. But the goal here is to show q-instability for almost all A, except for the case of A symmetric. If a non-symmetric matrix A is q-stable, then non-symmetry would imply that A_lj0 for some jl. For all such (j,l) pairs, two of the terms in equation 178 drop out. Dividing what remains by A_ljproduces:

[0509] (A_jl/A_lj)(RA^T)_ki=(RA^T)_ik  (179)

[0510] for all i and k. In order to satisfy this equation with A_jl=0, for all i and k, it would be necessary that (RA^T)_ik=0 for all i and k. This would require that R be singular (which, applied to equation 176, would violate the requirement that A be positive definite) or that A=0 (which is a symmetric case). Thus, A_jl0, and at least some combination of k and i exits such that (RA^T)_ki0. For every such combination of k and i, equation 179 may be divided to produce:

[0511] A_jl/A_lj=(RA^T)_ik/(RA^T)_ki  (180)

[0512] We may also deduce that the numerators and denominators of equation 180 are all nonzero, for these combinations (l,j) and (k,i). Because they must all be nonzero, we can apply the same logic for l and j interchanged, or for k and i interchanged, independently. This implies, for example, that:

[0513] A_jl/A_lj=A_lj/A_jl=±1  (181)

[0514] Furthermore, since we may use the exact same logic for any nonzero A_μν (with μ ν), using the same combination of k and i, we may deduce the following. For q-stability, for any matrix A which is not symmetric, every offdiagonal term would have to obey A_lj=−A_jl. In that case, equation 180 would also require that RA^Twould have to be completely antisymmetric.

[0515] This analysis could be extended further, to provide additional restrictions or requirements for q-stability. However, q-stability is possible here only in two situations—where A is completely symmetric, or where the offdiagonal terms are completely antisymmetric. Both situations represent sets of measure zero in the space of possible matrices A.

[0516] The findings above suggest the continuing use established Critic adaptation methods, and proving theorems which exploit the requirement that the schedule of experience x(t) should come from an appropriate probability distribution [73].

[0517] The Generalized Moving Target (GMT) problem previously discussed in [17] is discussed below. It will offer new gradient-based neural network designs to address that problem. It will also address related possibilities for assuring convergence more reliably (but still cheaply) for large-scale economic equilibrium models like the one discussed in [16]. However, it will also show how the application of these designs to the Critic adaptation problem here involve excessive cost, just like the previous stable methods of Landelius [43] and Bradtke [44].

[0518] The dynamic equations for M given above are all examples of a larger class of systems called Linear Differential Inclusions (LDI) [79]. In essence, the large problem to be considered is: for the system

[0519] ∂_tx=A(α)x,  (182)

[0520] starting from some bounded initial state x(0), prove that x(t) at later times is also bounded, regardless of the choice of α(t) at later times t. Equation 182 is not a linear dynamic system, because it refers to a family of matrices A(α), for different possible values of the parameter(s) α, rather than a single matrix A. The state vector being studied is actually the matrix M, and the set of parameters α is actually the state vector of the plant being controlled.

[0521] Boyd et al[79, p.74] cite many examples of systems of this sort which are unconditionally stable, but not q-stable. However, these examples give a complete intuitive feel for the nature of the gap between q-stability and unconditional stability.

[0522] Instead, suppose that x is simply a two-component vector, (x₁, x₂). Suppose that the parameter α takes on only two possible values. (α=1 and α=2, with: $\begin{array}{cc}A\left(1\right)=\left[\begin{array}{cc}a& 1\\ -1& b\end{array}\right];A\left(2\right)=\left[\begin{array}{cc}b& 1\\ -1& a\end{array}\right]& \left(183\right)\end{array}$

[0523] where a>0 and b<0 are both very small (near zero). In order to determine whether this system can be destabilized, how should α(x) be picked in order to maximize future |x|²over future times t. In this case, the choice between α=1 and α=2 is easy because:

[0524] ∂_t(½|x|²)=x^T(∂_tx)=ax₁²+bx₂²(for α=1) or bx₁²+ax₂²(for α=2)  (184)

[0525] Thus at each point x, more increase (or less decrease) in length is achieved by choosing A(1) when |x₁|>|x₂| and A(2) when |x₂|>|x₁|. If a is very small, and b is set such that b=−ka, and x is considered to be on the unit circle (without a loss of generality here), equation 184 is integrated over the interval from θ=−π/4 to θ=+π/4, in order to evaluate the total net increase or decrease in the length of x over that quarter circle. (The other quarter circles would involve the same increase or decrease, by the symmetry here.) In order to ensure a decrease of the length of x (i.e., in order to ensure stability), it is required that: $\begin{array}{cc}k\ge \left({\int }_0^{\pi /4}{\cos }^2\theta d\theta \right)/\left({\int }_0^{\pi /4}{\sin }^2\theta d\theta \right)=\left(\pi +2\right)/\left(\pi -2\right)& \left(185\right)\end{array}$

[0526] (Approximately k≧5 should ensure stability.)

[0527] It is easy to see that this system is not q-stable for any choice of k. By symmetry, we would expect Λ=|x|²to be the only serious possibility for a quadratic Liapunov function here. But at θ=±π/4, etc., both of the matrices A(1) and A(2) lead to increases in the length of x!

[0528] This example shows how a kind of negative boundary zone between A(1) and A(2) prevents the existence of a global quadratic Liapunov function. Nevertheless, if the dynamics of the system are guaranteed to move the system out of this boundary zone, then there may be no possibility of instability in the overall system. On the other hand, this example shows that a small destabilizing term like a (in A_jl(1)) requires a much larger counterterm (like −5a) elsewhere, in order to permit stability.

[0529] In a similar vein, for methods like HDP, consider the possibility of a dynamical system: $\begin{array}{cc}A\left(1\right)=\left[\begin{array}{cc}-1& 0\\ a& -1\end{array}\right]& \left(186\right)\end{array}$

[0530] where once again a is very small. Because A is not symmetric, HDP is not q-stable in adapting to this system. An algebraic analysis does show that there are training vectors like (1, ±½a) that can destabilize it, for some matrices M. However, continued training in such special cases would be expected to change M. In effect, there are negative boundary zones here as well, where M can grow, but only for a short time. Intuition suggests that HDP might well be unconditionally stable, in fact, for such a system A, even though it is not q-stable. Nevertheless, as with equation 184, these mechanisms are expected to work only for plants A which are close to q-stable plants (symmetric A) in some sense. This mechanism does not seem likely to prove that HDP is a universal stable adaptation rule.

[0531] Consideration of this example does give some clue as to why large interest rates may tend to stabilize HDP in some situations. Interest rates add a negative term to the diagonal of A, in effect. This cannot make a nonsymmetric matrix symmetric! However, it can make the matrix closer to a symmetric matrix, insofar as the nonsymmetric offdiagonal terms may become smaller relative to the diagonal terms.

[0532] The Generalized Moving Target ((GMT) problem may be defined as follows. Find the vectors v and w which simultaneously meet the two following conditions:

[0533] (1) E(v, w) is minimized as a function of v for w held fixed; and

[0534] (2) w=v.

[0535] In 1990 [17], it was shown that the equilibrium of HDP solves that problem, locally (meeting the first-order conditions for a minimum), for the error function: $\begin{array}{cc}E\left(\underset{\_}{v},\underset{\_}{w}\right)=\sum _t{\left(\hat{j}\left(\underset{\_}{R}\left(t\right),\underset{\_}{v}\right)-\left(\hat{j}\left(\underset{\_}{R}\left(t+1\right),\underset{\_}{w}\right)+U\right)\right)}^2& \left(187\right)\end{array}$

[0536] Special-purpose quasiNewton approaches may also be utilized.

[0537] The problem with stability, using HDP, may be explained as the result of minimizing E with respect to v, while not anticipating the effect that this has on changing the problem by changing w. But minimizing E(v,v) directly, as in the Galerkinized methods, would solve the wrong problem. Likewise, all manner of Bayesian updates rooted in the effort to minimize E(v,v) would be nothing but sophisticated ways to solve the wrong problem, to move to an incorrect equilibrium.

[0538] One approach is to try to solve the GMT problem in general, and then apply the result to the adaptation of Critics. As a preliminary to finding a general solution, consider the quadratic special case of GMT, QGMT, where:

[0539] E(v,w)=v^TAv+w^TBw+v^TCw+x^Tv+v^Tw,  (188)

[0540] where A, B, C, x and v are arrays of parameters held fixed, The first-order condition for a minimum of equation 188 with respect to each component v_iof v is: $\begin{array}{cc}E\left(\underset{\_}{v},\underset{\_}{w}\right)=\sum _{j}2A_{\mathrm{ij}}{v}_j+\sum _j{C}_{\mathrm{ij}}{w}_j+x_j=0& \left(189\right)\end{array}$