Background of the Invention

The present invention relates to industrial control systems.

Industrial control systems are used to control a wide variety of industrial processes. Examples of such processes include the control of steel rolling mills (to maintain uniform output thickness of rolled steel strip ⁾ , ship stabili¬ zation, dynamic ship positioning, adaptive ship autopilots, temperature control of buildings, radio aerial tracking, and a wide variety of chemical processes.

In all such processes, the problem in general terms is to maintain a desired quantity (such as the heading of a ship) constant (or following a desired pattern) in the face of peculiarities of the apparatus itself (such as slow response time of the rudder in turning to a desired position, and of the ship in turning in response to a rudder position change) and external influences (such as steady cross currents and varying wind forces). The control system (controller, actuators, and sensors) measures the desired quantity and generates a control signal which is fed into the system being controlled (as a control signal to the rudder position gear).

The controller can be regarded, in mathematical terms, as calculating the control signal as some function of the measured quantity. The particular function chosen can generally be regarded as a transfer function having a parti¬ cular form ⁽ e.g. a polynomial of given degree) and having particular parameters ⁽ which would then be the coefficients of the various terms of the polynomial). The form of the function is chosen from a knowledge of the general character¬ istics of the system to be controlled. The parameters of the chosen form can then be determined in various ways, involving measurement or calculation of the characteristics of the system to be controlled.

In many control systems, the controller parameters are fully defined by the initial design, and the controller provides satisfactory control. However, in some systems, the operating conditions and/or system characteristics vary widely with time; for example, the characteristics of a cargo ship will differ widely depending on its speed, whether it is fully laden or in ballast, etc. In βuch

circumstances, it may be desirable for the control system to be adaptive. For that, the characteristics of the system are monitored and the parameters of the controller are adjusted accordingly. Thus the control system does not merely calculate the control signal in dependence on the measured quantity; it also adapts itself to changes which are detected in the general characteristics or behaviour of the system being controlled by continual parameter adjustment.

One known controller design technique is known as H» (H-infinity) design. This has various advantages, but it has the drawback that the calculations involved in determining an H® controller in any particular instance are in general very complicated and involved. The H» method has therefore generally been used for the initial (off-line) design of controllers.

The main purpose of the present invention is to provide a control system which is of the H» type (in a sense discussed in detail later and which is sufficiently simple in the calculations involved for the controller calculations to be carried out in an on-line system. As explained above, such a system is an adaptive system, because the controller is being updated regularly.

Accordingly the present invention provides an industrial controller for controlling an industrial process, comprising a controller unit fed with a reference signal (r> and the output (y> of the process, a parameter estimator fed with the output <u> of the controller unit and the output (/> of the process and producing estimated polynomials A, B, and C, and a controller design unit fed with A, B, and C and feeding the controller unit with coefficients for implementing a polynomial ratio control function C _n / ), characterized in that the controller design unit comprises means for storing two functions F _e <= F _cr , F _ceJ ) and P _c (= P _e -,/P _ee _ ⁾ , means for calculating the unstable zeros L_ of the function L (= Per. Fees B - F _e „ P _eeJ A ⁾ , means for solving the matrix equations

F A P _c «_ λ + L_ G = P _en C F* r"

F B F _ed λ - L H = F _en C F* zr", and means for calculating Cn and C _d as G F _ce. and H P _cei .

The equations may be solved either by an eigenvalue/eigenvector calculation or by F-iteration on the equations

F A ? _ees λ + L_ G = P _cr , C

F B F _ed λ - H = F _cr , C

in which F* is calculated as the adjoint of F and the right-hand side of these equations are multiplied by F*.

The parameters of the stored functions P _c and F _c may be adjustable.

Listing of Drawings

Two embodiments of the present invention will now be described in detail, by way of example, with reference to the accompanying drawings in which:

Fig. 1 is an idealized diagram of a ship;

Fig. 2 is an abstract block diagram of the heading control system of the ship;

Fig. 3 is a more detailed diagram of a control system embodying the invention;

Fig. 4 is an idealized diagram of a steel rolling mill;

Fig. 5 is an idealized gain/frequency graph for the rolling mill;

Fig. 6 is an idealized gain/frequency graph illustrating a cost function; and

Fig. 7 is a more detailed block diagram of parts of the Fig. 3 system.

Introductory Summary

The present invention provides an industrial control system which comprises a controller controlling an industrial process (such as a ship, a rolling mill, etc ⁾ , the controlled variable of which is y. The controller and process outputs are fed to a parameter estimator which estimates the operating polynomials A, B, and C of the process. The operating polynomials are transfer functions (of the Laplace operator ε if the control system is a continuous one, or of the unit delay operator zr~ if the control system is a sampled one. A controller design unit processes the output of estimator to produce coefficients for the con¬ troller, which implements a polynomial ratio C _n /C _d .

The controller design unit implements an on-line H» algorithm which is based on minimizing a "square of sum" form of function (P _e .u(t) + F _c .e(t)) ² .

(This results in equations which are easier to solve than those arising from the usual "sum of squares'* form of function. ⁾ A polynomial

L *■**• P _cr , F _cd B — Fen P _c <3 A is defined, and the unstable factor L_ is calculated.

The following polynomial equations arise from the optimization problem: F A P _ed λ + L G = P _cr , C F* z-" F B F _ed λ - H = F _er> C F* r-\

These equations may be solved by the eigenvector/eigenvalue technique. For this, the eigenvalue with the largest magnitude is extracted, the correspond¬ ing eigenvector is extracted, and C„ and C _β calculated as G F _c< -, and H P _cc , for passing to the controller unit.

Instead of this, an iterative technique (the F-iteration technique) can be used, in which the equations

F B F _cd λ - L H = F _er , C are solved, F* (the adjoint of F) is calculated, the polynomials on the right-hand side of these equations are multiplied by this adjoint, and the procedure is iterated.

F _c and P _c are chosen so that P _e dominates at low frequencies and F _c at high frequencies. The parameters of F _c and P _c may be adjustable, and P _c may be chosen to include a term which reduces the responsiveness of the system to disturbances in a frequency band in which noise is strong.

Detailed Deεription

Fig. 1 shows a ship 1 in idealized form. The ship's heading is measured as an angle ψ (with respect to a fixed direction, such as true North), and the position of its rudder 2 is measured as an angle δ (with respect to the ship's long axis). The ship is subjected to various disturbing forces, such as a wave action ⁽ where the waves may be of a certain frequency and amplitude and travell¬ ing generally in a direction 3 - all of course known only statistically) and wind action ⁽ where the wind direction 4, its strength, and its gustiness are also known only statistically). The object is to maintain a given heading, and this is to be achieved by a suitable control system.

Fig. 2 shows the steering control system of the ship 1 in abstract form. The ship and its steering gear are indicated by a box 10, the desired heading is given as a signal REF, the signal fed to the ship and steering gear is shown as a signal u, and the actual heading is shown as a value y. Ignoring the blocks 13 and 14 for the moment, the actual ship heading y (taken for the moment as equal to the output signal m from the block 10 ⁾ is measured and fed to a summing unit 12 which also receives the reference signal REF and produces the signal e which represents the "error" or difference between the desired and actual headings. This signal e is fed to the controller, represented by block 11, which produces a control signal u which is fed to the ship ⁽ block 10 ⁾ . (This particular representation of the process in terms of A, B, and C is known as the ARMAX model. ⁾

The feedback loop through the summer 12 and controller 11 is negative, as indicated by the minus sign at 12 for the feedback signal path, in the sense that any small disturbance of the system is cancelled out so that the desired heading y is held constant. One form of "disturbance" is a change of the reference signal REF. It is desirable for the control system to cause the output y to change to follow the signal REF, and to cause this change to occur as fast as feasible, but to avoid excessively large swings and oscillation of the output y in the course of the change.

The entire system can be represented or modelled in various well known ways, of which the way just described in one. In many instances, the behaviour of the system being controlled can be represented by means of linear differential equations. Thus for the ship shown in Fig. 1, this approach would involve developing a set of differential equations which describe the operation of the ship. This process would start with an equation relating the angular accelera¬ tion ⁽ rate of change of angular velocity) of the ship, its angular velocity (rate of change of yaw ⁾ , and the torque (dependent on the rudder angle and other hydrody namic factors), thus: τψ + ψ = Aδ. Further equations relating to the effects of wind and wave and various parame¬ ters would also be developed. (One important parameter is the ship velocity, affecting the values of k and τ, since the turning effect of the rudder depends strongly on the ship velocity.)

It is conventional to take the Laplace transform of such equations; this yields what are known as transfer functions. Thus each of the two blocks 10 and 11 can be represented as a transfer function, and for present purposes a transfer function can be taken as the ratio of two polynomials- The transfer function for the block 10 is shown as A/B, where A and B are polynomials, and that for the block 11 is C^/Cc, where C-, and Ca are polynomials forming the numerator and denominator respectively of the transfer function of block 11. (It should be noted that the C„ and C* of block 11 ^" are distinct from the C of block 13. C-, and Co are computed by a controller design module in the controller, while C (or C/A) is determined by the nature of the system to be controlled. ⁾

If the system is more naturally regarded as a sampled or discrete-time system, then its representation is by difference equations instead of differential equations. In such cases, it is conventional to take the z-transform of the equations; this is somewhat analogous to the Laplace transform, in that it yields functions which behave analogously to the transfer functions of continuous systems mentioned above and are described by the same term, transfer functions.

A complication which often occurs in practice is that the system may include sources of disturbance having characteristics which are not fully known. In the present example, the effect of crosswinds on the ship form such a source of disturbance. Such a disturbance can be represented as a noise source. Often the noise is not white, but is e.g. strongest in some particular frequency band; it is conventionally represented by a white noise signal w which is passed through a colouring filter 13 ⁽ which is a transfer function which can be chosen in the form C/A, to match the denominator of the transfer function of block 10) producing the actual wind noise signal d which is combined with the output m of the block 10 by a summer 14 to produce the actual heading signal y. Further disturbances can similarly be included by including further summers in series in the feedback path from the output of the ship 10 to the refernce signal summer 12.

The general objective of the control system is, as noted above, to achieve "good" control of the system, in the light of the system characteristics (trans¬ fer function B/A ⁾ and the noise signals. This involves computing suitable poly¬ nomials C-, and C _α .

One way of solving this problem involves a once-for-all calculation of the feedback transfer function. If the system being controlled has well-understood and relatively simple behaviour and noise effects are small, its transfer function can be calculated (using calculated moments of inertia, for example), and the optimal controller calculated. If the system being controlled is not understood too well or cannot have its transfer function calculated readily (for example, because the moment of inertia of a part with a complicated shape is required), then it may still be feasible to operate the system for an initial period during which its characteristics are measured, and an optimal controller for the system can then be calculated on the basis of those measurements.

If the system being controlled has a significant input noise disturbance, then the control transfer function is preferably chosen to minimize its effects, e.g. by producing a strong corrective action at the frequency band where the noise is concentrated (so as to compensate for the noise) but weaker action at other frequencies. The wind disturbance discussed above is input noise, i.e. noise arising "inside" the process being controlled; it is shown as entering the signal chain after the process 10 merely for convenience of analysis.

The system may also be subject to an output noise disturbance (not shown), which would act between the summer 14 and the controller block 11. Such a disturbance could arise for example from measurement errors in measuring the actual heading y, and would result in the signal fed to block 11 being different from the actual heading. The controller should be designed to neglect such output noise disturbance.

_ If the characteristics of the various noise signals are substantially constant, the once-for-all off-line design procedure may still be satisfactory.

The calculation of the optimal controller usually involves the minimization of a function, known as the "cost-function", which is a function of the system characteristics ⁽ including the controller transfer function) and a function defining the various features which are desired to be minimized (e.g. rapid response, low overshoot, good stability, etc) and their relative weights. The cost-function is discussed in more detail later.

One approach to optimal control is known as the H» technique, and involves a well known procedure of minimizing the εupremum or maximum value of a

frequency response which can be defined by system sensitivity and complementary sensitivity functions. The numerical computation of H» controllers is generally very difficult. Such controllers are therefore normally computed only by an off-line procedure.

In more complicated systems, the above off-line approach is not satisfac¬ tory. For example, in the case of the ship heading control system discussed above, the wind noise characteristics (amplitude and frequency spectrum) may vary considerably as a function of time, and similarly the ship characteristics will vary depending on its loading and speed. In such cases, some form of adaptive or self-tuning system is desirable. In such a system, the system parameters are calculated repeatedly on the basis of the immediate past history of the system (i.e. the values of the input and output signals of the system over a suitable period up to the current time). This process is known as system identification or parameter estimation. The controller transfer function is updated regularly on the basis of the calculated system parameters, and used to control the system.

It is important in such on-line adaptive systems for the parameter estima¬ tion and controller calculation algorithms to be reasonably simple and straight¬ forward, with guaranteed convergence properties, in order that the control calculation is simple, reliable, and completely automatic.

A major problem with many self-tuning or adaptive techniques is that they are based on the assumption that the system being controlled can be represented by a relatively simple model, - i.e. a model in which the two polynomials A and B are of relatively low order. This assumption is often false. Thus although the self-tuning system adapts to variations in the low-order approximation to the actual system being controlled, it does not take into account the fact that the model is only an approximation to the real plant.

One way of describing a system which is to be controlled is by means of a graph plotting the gain of the system against frequency. More specifically, the gain is shown as the Y co-ordinate against the frequency, plotted linearly along the X axis. The "gain" is the ratio (usually expressed in decibels) of the response of the system (measured as some suitable output quantity) to an input driving signal. ⁽ This form of description is commonplace for electrical systems such as amplifiers. ⁾ This type of description does not encompass all the prop-

erties of the system; for example, it fails to describe phase εhifts ⁽ which are often shown as further curves on the same graphs.) However, gain graphs give a good general indication of the nature of the system.

Fig. 4 is a typical gain/frequency graph for an idealized steel rolling mill shown in Fig. 5. The mill consists of two rollers 40 and 41, between which steel strip is passed to be rolled to a desired thinness. Variations in such variables as the input thickness, temperature, and consistency of the eteel strip and the speed of the rollers affect the output thickness, and the speed of the rollers can be adjusted accordingly. The rollers are driven by means of two motors 42 and 43, which are connected to the rollers 40 and 41 through shafts 44 and 45.

The response of the system to variations in the input (control) signal is strongly dependent on the frequency of the input signal change. This is because the system incorporates a number of inertias, such as the inertia of the rollers 40 and 41 and the inertia of the steel strip being rolled. The gain/ frequency graph therefore generally slopes downwards to the right, as shown. The curve has various points at which its slope changes, such as points 50 and 51; these points correspond to corresponding linear factors in the transfer function polynomials A and B.

The low frequency response can be modelled reasonably well, from a know¬ ledge of the mechanics of the system (i.e. a knowledge of the various inertias, etc.). Such modelling will in general not be perfect, because some properties of the system will not be known accurately; for example, the inertia and consistency of the steel strip being rolled will vary from strip to strip. These uncertain¬ ties can be described as parameter uncertainties, since the details of the model will depend on such parameters as the inertia and consistency of the steel strip being rolled. Such parameter uncertainties can be dealt with by using an adapt¬ ive control system.

Fig. 4 also shows some high frequency variations 52 in the response graph. These variations are due to such matters as resonances of the roller-40/ shaft-44/motor-42 and roller— 41 /shaft-45/motor— 43 combinations. (The motors 42 and 43 are normally housed in an enclosure at some substantial distance from the rollers 40 and 41, so the shafts 44 and 45 are of substantial length, and such resonances are therefore substantial.) Although the occurrence of such high

frequency features can be expected, it is often impracticable to attempt to model them, because their precise nature may not be calculable in practice; they may be termed unmodelled dynamics or unstructured uncertainties.

Such high frequency features can often be ignored; as a general rule, they will not rise up as far. as the 0 dB gain level, and ignoring effects below that level will not have disastrous effects. However, it is possible for such effects to approach the 0 dB level. Hence it is desirable to have a controller which is robust in the presence of such uncertainties. It is not feasible to do this by adaptive techniques. Instead, such compensation can be achieved by H« tech¬ niques. It is therefore often advantageous to apply both adaptive techniques and the H» technique to a given single control system, since the two techniques affect separate aspects of the system: adaptive techniques perform parameter adjustment of low-frequency modelled features while H« techniques compensate for unmodelled high-frequency features.

The H» technique mentioned above allows for uncertainties in system des¬ cription. However, this technique has hitherto been used only for designing fixed (off-line) controllers. One important reason for this is that, as noted above, the calculations involved in using the H» technique are very complicated, and cannot conveniently be utilized in an adaptive or self-tuning system.

In the present invention, the two design approaches are combined. This results in the best features of a self-tuning system being combined with those of an H« design. Thus the present invention provides a self-tuning controller which identifies a low-order approximation to the system being controlled but which takes into account the modelling uncertainties in the approximation. The resulting system combines the adaptiveness of self-tuning with the robustness and frequency response shaping provided by the H« technique.

Considering robustness in more detail, it is known that this can be improved in fixed (non-adaptive) controllers by an appropriate choice of frequency depen¬ dent cost function weightings. This characteristic can be maintained in the present invention. For example, good sensitivity characteristics at low fre¬ quency can be achieved by introducing an integrator in the cost function on the error term. This ensures that the controller has high gain at low frequency, to ensure that the system tracks low frequency reference signal changes closely and also minimizes the effect of low frequency input disturbances. The cost

function weightings can also be chosen to minimize the effect of high frequency dynamics which are not adequately represented in the model, and which are a source of instability in known self-tuning controllers. Known self-tuning controllers simply minimize a least-squares cost function; while this is good for disturbance rejection, it has adverse effects on robustness.

A further problem with existing self-tuning controllers is that most of them have poor performance with non-minimum phase systems. A minimum phase system is one in which the transfer function polynomials have no factors with positive real parts, and the model normally involves this assumption. But many actual systems are not minimum phase systems; in particular, systems involving sampling often behave as non-minimum phase systems. This results in poor per¬ formance of existing self-tuning controllers when applied to such systems.

The techniques of the present invention can also in principle be enlarged for very complex systems, including multi-loop systems and multi-variable . systems.

Fig. 3 shows a control system, such as the system of Fig. 2* in a more detailed form. Block 20 represents the process being controlled, and corres¬ ponds to block 10; unit 21 represents the controller, and corresponds to block 11 together with the summer 12 of Fig. 2. (The controller 11 of Fig. < '_ can be located in the path between the summer 12 and the process block 10 instead of where shown. Such a change of location is not of great significance; there are well-known techniques for converting the transfer function for the controller from the form suitable for either location to a form suitable for the other, and the two forms have generally equivalent characteristics.) Noise sources are not shown in Fig.3, although they will of course normally be present.

Block 20, the process, is fed from unit 21, the controller, via a digital-to- analog converter 30 and an actuator 31, which controls the physical variable of the process - e.g. the position of a rudder. The output of block 20 is another physical variable - e.g. a ship heading - and this is measured by means of a sensor 32 and converted to digital form by an analog-to-digital converter 33.

The controller unit 21 performs a simple control function which applies the transfer function C--/C-J to the difference e between the feedback signal y and the demand signal r ⁽ signal REF of Fig. 2). This transfer function is not fixed

but is constantly updated. The control iε achieved by a two-etage process. The control signal u from the controller 21 to the process 22 and the process output variable signal y from the process 20 are fed to a parameter estimation unit 34, and the parameters estimated in unit 34 are fed to a controller design unit 35, which uses them to "design" the controller 21 - i.e. continuously update the two polynomials C„ and C _d of the controller unit 21.

The estimation unit 34 utilizes known principles, as does the controller unit 21. It is known to provide, between these two units, a controller design unit which incorporates a controller calculation algorithm. Examples of such algorithms are LQG (Linear Quadratic Gaussian) algorithms, minimum variance algorithms, etc. (An H» controller is somewhat similar to an LQG one, but the LQG controller minimizes the average deviation of the error over frequency, while an H» one minimizes the maximum deviation of error over frequency. The H» closed loop system iε therefore more robust, i.e. less liable to result in osc¬ illations.) It is in unit 35 of the present system that the H» type calcula¬ tions are performed.

The estimator unit 34 receives, at each sample instant, the value of u cal¬ culated on the previous sample instant (or possibly two or more sample instants previously, depending on the delay in the plant) and the current measured value of y. Using the data received at each sample instant, the estimation unit uses a recursive method to build up a transfer function representation of the system - that is, to identify the plant parameters or polynomials A, B, and C. The estimation unit may be preset with initial values for these functions to use on start-up, to minimize the time required for the estimated values to converge to a reasonable approximation of the actual values.

The estimator has just been described in terms of a sampling system, and the following description is alεo primarily in terms of a sampling system and the unit delay operator z -Cor z~ ^' As noted above, however, the same principles can be utilized with little change in a continuously-operating system and using the Laplace operator ε.

The fundamental purpose of the control system is to cause the process var¬ iable y to follow the demand signal r. The tracking error e ⁽ t ⁾ is the differ¬ ence between these two signals, and the minimization of this tracking error iε therefore a fundamental desideratum. However, it iε also desirable to minimize the control action u(t> which iε required to achieve the tracking of the process being controlled to follow the demand signal. This iε because if the control action is large, stability can be adversely affected and the physical demands on the actuator 31 and the plant in which the process is being carried out can become excessive. There is thus a conflict between these two desiderata. To achieve a satisfactory balance between these two conflicting desiderata, a function to which both signals contribute iε therefore selected and minimized. This function is normally defined, in H» controllers, by using a first weighting function P _c for the weighting of the error signal and a second weighting function F _e for weighting the demand signal.

The operation and function of substantially any controller can be described in mathematical terminology. In particular, the operation and function of the present controller can be described in terms of the equations

F A P _cd λ + L_ G = P _cr , C F* z-"

F B F _ed λ - L H = F _en C P r"

The present controller can be regarded as determining the minimal degree solution of these equations.

In these equations, P _cn and P _ce , are the numerator and denominator of a function P _c , F _c „ and F _ct _ are the numerator and denominator of a function F _c , and F and F* are adjoint functions (i.e. F (z ^-'1 ) = F* (z)>. n (or more explicitly n _r ) is the degree of F; the polynomialε on the left-hand sides are in powers of z ^{- •} , and the multiplicative factor z~" iε required to ensure that the right-hand sides are also in powers of z~ ^~ rather than of z. L_ is obtained from a polynomial L which is defined by

L - P _{en cc} ι B - F _cr> P _ee) A.

This polynomial is of course factorizable into a eet of linear factors of the form z - α, where the α'ε can be real or complex (the complex roots (zeros) occur, of course, in conjugate pairs). The α's are the zeros of the polynomial,

and are termed unstable or stable according as their moduli are or are not greater than 1. (That is, an α is an unstable zero if lαi > 1. ⁾ L_ is the product of all unstable factors of L, i.e. the factors with unstable zeros.

A, B, and C are of course the polynomials discussed above, The manner in which the functionε P _c and F _c are choeen is discuεεed below, It should be noted that F, is distinct from F.

The preeent controller can be described in mathematical terms as: (a) solving these equations by determining the unknowns λ (which iε a scalar) and F, G, and H (which are polynomials of minimum degree in z), and; (b) from this solution, determining the transfer function C _n /Ce - (G F _cej ) / (H P _ct3 ) . This transfer fuction iε implemented by the present εystem.

The■ plant iε assumed to be linear. The conditions of this assumed linear plant may change, either because the plant is in fact non-linear and its operat¬ ing point changes or because there are changes in its load or other operating conditions. The model estimated will change correspondingly, and the controller will be recalculated accordingly. Thus the adaptation takes place continuously, and so the present system provides a self-tuning controller.

In the present εystem, there are two different ways in which these equa¬ tions may be solved. One is the eigenvector-eigenvalue method, and the other is the F-iteration method.

In the eigenvector-eigenvalue method, the two equations are written in matrix form

f CP _e „ C) -tCL) 0 F t(A P _ce _> 0 0 F

- λ t' <F _β „ C) 0 -t(L) G tCB F _ee _> 0 0 G = 0 H H

Each of the polynomials F, G, and H is converted into vector form by having its coefficients written in a column, and the three columns are written one above the other in the second matrix of each product in the equation. The polynomials Pe . C etc ⁽ which are, of course, each a product of two polynomials) are converted into matrix form b bein written as Toe litz matrices.

The Toeplitz matrix t (X) of a polynomial X iε a matrix with the first column formed by writing the coefficients of X downwards (starting with the lowest order (constant) coefficient), and each successive column εhifted downwards one step relative to the previous column, the number of columns being determined by the number of rows in the matrix which the Toeplitz matrix iε being multiplied by; spaces above (for all but the first column) and below (for all but the last column) the coefficients are filled with 0*s. The Toeplitz matrix t' CO is a matrix like t (X) but reflected in a vertical axis (for the discrete caεe, uεing the unit delay operator z) or with every element oc in all even-numbered columns being replaced by its complement -α (for the continuous caεe, uεing the Laplace operator ε). That iε, the Toeplitz matrices are

a ₀ 0 0 0 a- a ₀ 0 a ₂ a, a ₀ t (A) =

0 0 0 0 0

or the three matrices respectively.

This matrix equation is in the form

(Q - λR> S = 0

where Q and R are known matriceε, S iε an unknown vector, and λ iε an unknown scalar. This matrix equation is known as the generalized βigen problem, which is a well-known problem in linear algebra. The matrices Q and R are square matrixes of size n * n, and S iε a column matrix corresponding to the polyno¬ mials F, G, and H. This equation yields n solutions for λ and S. In the present εystem, it iε neceεεary to choose the solution which yields a stable control system. This iε achieved by choosing a solution with an extreme value of λ - either the minimum or the maximum, depending on the formulation of the eigen problem. Once the equation has been solved and the polynomials F, G, and H found, the transfer function C _r , C _d iε calculated as noted above, and this transfer function is passed from unit 35 to unit 21.

Any method of solving the eigen problem can be used. Certain known methods of solving the eigen problem, such as the QZ method, have the charac¬ teristic of operating in two stages, the first stage being the determination of .the eigenvalues and the εecond being the determination of the corresponding eigenvectors. In uεing such a method, the eigenvalues would first be calcu¬ lated, the extreme one chosen, and then the corresponding eigenvector calculated.

Such a method allows a refinement, since the method allows the calculation of eigenvalues only. For this refinement, the εyεtem monitorε the way the value of the chosen eigenvalue changes, and calculates the corresponding eigenvector only if the eigenvalue changes significantly. This reduces the average calcula¬ tion load, enabling the εyεtem to follow the eigen value more closely during periods when it is changing only slowly. Alternatively, the system can be used for unrelated background processing during such periods.

In the F-iteration method, the H» two equations above are solved by solving a series of similar but simpler equatione. The simpler equations are F A P _ces λ + L_ G •= P _cn C F B F _ed λ - L H = F _cr , C

These are in fact the LQG equations, which represent the LQG controller for the εame εystem; as with the H» problem, the LQG controller iε given by the transfer function ⁽ H P _ed ⁾ -' ⁽ G F _c< -, ⁾ . In matrix form, these equations can be written as

This equation ie in the form Q R = S, which ie a εtraightforward syεtem of linear equations, in which Q and S are real matrices and R iε a real vector, and the unknowns are the elements of R. This system of equations may be readily solved by any of a variety of well-known techniques. When R has been found, i.e. the polynomials F, G, and H are found, then F*, the adjoint of F, is calculated. The polynomials on the right-hand side of these equationε are multiplied by this adjoint - i.e. the right-hand sides of the equationε are replaced by the products P _cτ - C F* and F _cr C F* . The equationε are re-εolved, and the process is iterated. It can be shown that the iteration converges to a solution of the equations given above.

In practice, in a sampling, system the number of iterations performed during each sampling period may be any number from 1 upwards. If the number of iter¬ ations iε large enough, then a "pure" H» controller of the above type iε imple¬ mented. If the number of iterations is small, then the controller is a mixed controller with LQG characteriεticε as well as the characteristics of the above type, but is always converging towardε a pure controller of the above type;** its LQG content will be the result of only recent changes in the plant being con¬ trolled, and will always decay to nothing if the plant conditions remain stable.

Fig. 7 is a more detailed block diagram of the present controller, showing the details of blocks 21 and 35 of Fig. 3; The components of a sampled or discrete εystem are shown.

Considering first the controller, block 21, this implements the equation

u = C _π /C _d X e

C„ and C _d are polynomials in r ¹ , say

C„ = q ₀ + q..z~ ^τ + q-≥.z- ² ,

Inserting these into the original equation and rearranging it to obtain u gives

p ₀ . u = t Cq ₀ + q. .z~' + t _z .__- ^~ ~').e_/i . * _i .z- ^~ + τ..__-~ ). ιΔ.

This gives the present value of u in terms of the present and past values of e and the past values of a This is implemented by means of the circuitry 60-66. The signal e iε fed to a chain of unit delay circuits 60, whose outputs are fed through a eeries of scaling circuits 61 which multiply the various e values by the q coefficients to an adder circuit 62. Similarly, the signal u is fed to a chain of unit delay circuits 63, whose outputs are fed through a series of scaling circuits 61 which multiply the various u values by the p coefficients to an adder circuit 65. The outputs of the adders 62 and 65 are fed to a division circuit 66, which produces the signal u by dividing the output of adder 62 by the output of adder 66.

It will be noted that there is no scaling circuit for the undelayed υ signal, as this undelayed εignal does not appear on the right-hand side of the equation for u. Further, it iε assumed that the coefficients have been scaled so that po *= 1; this can readily be achieved in the calculation of these coeffi¬ cients. The scaling circuits 61 and 64 are adjustable or programmable to oper¬ ate with varying scaling factors; as will be seen, the p and q coefficients are supplied by block 35 and are variable.

The function of the parameter estimation unit 34 is to generate the poly¬ nomials A, B, and C which were discussed with reference to Fig. 2. This estimation is performed in a known manner, as described for example in L Ljung and T Sδderstr ^β m, "Theory and Practice of Recurεive Identification", MIT Preεε, London, 1983. More precisely, the ⁽ maximum) orders of these polynomials are chosen by the designer, and unit 34 generates estimated values for their coeffi-

cientε. Thiε is performed recursively, uεing a set of equations which update the current estimated set of polynomials, i.e. sets of coefficients a,,,, b _m , and c,„, by adding a set of correcting values which are calculated to minimize a predic¬ tion error. The prediction error is closely related to the "residual", which iε the value of the noiεe signal as determined by the estimated polynomials, the known past values of the signals u and y, and the estimated past values of the signal w (Fig. 2). Estimated values of the noise signal have to be used because the actual values cannot be measured directly; the values have to be obtained by calculation, and the calculation involves uεing previouεly estimated values. The estimation algorithm includeε a "memory length" or "forgetting factor" parameter which determines how many past values of the variables are required by the algorithm.

The controller design unit 35 receives the estimated polynomials A, B, and C (that is, the coefficients of these polynomials) from unit 34. It also has built into itself a pair of weighting functions P _c and F _c (each in the form of numera¬ tor and denominator polynomials). From these, it calculates the polynomials C-, and C , i.e. the coefficients of these functions (the p's and q's), and passes these coefficients to the controller unit 21, where they are uεed to εet the weights of the weighting units 61 and 64.

In more detail, the unit 35 has a store 70 containing the polynomials P _cτ -, Pc _d t F _cn , and F _cel , and a set of input ports 71 at which the polynomials A, B, and C appear. (All polynomials are in the form of sets of coefficients.) The store 70 and ports 71 feed a multiplier 72 which calculates the polynomial L, in accordance with the equations above. This feeds a factorizing and dividing unit 73, in which the factorε of L are determined in a well-known manner (the numer¬ ical extraction of factorε from polynomials is conventional). The factorε (which are in general complex) are classified according to whether or not their moduli are greater than 1, and unit 73 divides out from L those factors with moduli less than 1, leaving a reduced polynomial L_ including all factors with zeros whose moduli are greater than 1. (An analogous situation ariseε in the contin¬ uous-time case.)

The polynomials from memory 70, ports 71, and unit 73 are fed to a matrix unit 74, which combines them into the two matrices in the matrix equation above. This unit is coupled to an eigenvector extraction unit 75, which extracts the eigenvectors λ and determines the largest eigenvalue. The largest eigenvalue is

paεεed to an eigenvector unit 76 which determinee the correεponding eigenvector from the matrix equation held in unit 74. This eigenvector incorporateε the two vectorε G and H, and theεe are passed to a multiplier unit 77 where they are multiplied by the polynomials F _ceJ and ? _et _ to produce the polynomials C,, and C«a, which are passed to the controller unit 21.

If the F-iteration method is used, then the matrix unit 74 generates the simpler matrix equation discussed above, and the eigenvalue unit 75 and eigen ^¬ vector unit 76 are replaced by a linear equation unit which generates the poly¬ nomial F*, the adjoint function discuεεed above, which iε fed back to the matrix unit 74 to produce an iteration of the F equations, as discussed above, and iter¬ ation proceeds as also discuεεed above. The linear equation unit also produces the polynomials G and H at each iteration, and after the iteration haε been completed, these polynomials are fed to unit 77 as before.

The variouε calculations and manipulations of information may, be performed by meanε of diεtinct circuit elementε aε illuεtrated, or to any convenient extent by meanε of computer means which is suitably programmed to εimulate the blocks and functions shown.

Depending on the particular application, there may be several iterations per sampling period, a single iteration per sampling period, or even several εampling periods per iteration. The number of iterations per sampling period will normally be fixed. However, there may be circumstances in which the number is variable. For example, if the conditions of the syεtem normally change only slowly but there are easily recognized situations in which they are known to change rapidly, it may be appropriate for the number of iterationε per sampling period to be increased temporarily in such conditions.

It may be noted that depending on how elaborate the model of the εyεtem iε, the number of equationε iε not necessarily two. For example, if a more general model with more than one noiεe input is assumed, there may be three equations.

The present invention is a form of H« controller, which is defined in the following way. The basic principle of the H» controller iε that an optimal control criterion iε minimized and iε expreεεed in termε of a cost function J defined as:

J = su ⁽ X )

where sup over lzl=l repreeentε the mathematical εupremum over all frequen¬ cies, and iε the norm associated with the H* normed space of complex functions. The function X is a cost index, the cost function of which is to be minimized.

The function X depends on the feedback controller. There is therefore eome controller which minimizeε X εuch that X haε a constant value over all frequencies. In this case, the controller is an H» controller and is unique, since there is only one controller which minimizes the cost function.

It can be shown that the present εyεtem representε a minimization of a function X which iε the power εpectrum of a weighted sum of the tracking error e (= r - v> and the controller output εignal u, as follows:

X = power εpectrum of ( P _c (z ^-1 ) e(t> + F _e (z ^-* ) u ⁽ t ⁾ >

where P _e and F _e are defined in the eame way as in the standard H» εyεtem. Thiε function iε a "εquare of sum" function, in contrast to the more usual "sum of squares" function employed in both LQG (Linear Quadratic Gaussian) εyεtems and standard H» εyεte ε.

In more detail, when an optimal control εtrategy iε being deviεed, εome mathematical criterion which represents a "cost" iε choεen, and this "cost" must then be minimized. For example, if we are concerned with a motor vehicle, we might choose the "cost" aε the quantity of fuel uεed in a given period of time. This will be high for very low εpeede (becauεe of the low gearing and high engine epeed) and for very high εpeeds (because of engine losses at high revs and wind resistance). So there will be some speed at which the fuel used in the given time iε a minimum; thiε is then the optimal speed, and the correspon¬ ding fuel consumption is the optimal fuel consumption. The principle of .the optimal feedback controller calculation iε then analogous to taking any car and, from its specification, calculating the speed at which the petrol conεumption iε a minimum.

To return to the optimal feedback control syεtem described earlier, it iε convenient to chooεe, aε a cost-function, a function of the tracking error e and the control εignal u It is desirable to minimize both these quantities. Minimization of the tracking error is desirable for obvious reasons, and minimi-

zation of the control εignal iε desirable so aε to minimize the size of actuatorε required and their power conεumption.

The general form of the cost-function is then:

T J = J F (e(t>, u(t) > dt 0

where F iε the coεt index, and e and u are dependent on the controller. There will be one (or poεεibly more than one) controller which minimizes J; this iε the optimal controller. F (which iε here a function, distinct from both the poly¬ nomial F and the function F _e ) is selected such that a desired criterion which involves e and u iε minimized.

The principal form of the coεt-f unction iε a εum of squares function, in which the cost index being minimized iε a εum of squares of the variableε - in thiε caεe e and υ. Thus thiε cost-function is typically

^• T J, = J (P _c . e ^z + F _c . u ² ) dt 0

Thiε iε analogouε to the familiar use of a mean εquare function in such matters aε obtaining a measure of the magnitude of an oεcillatory εignal (where the εquare root of the mean εquare iε in fact uεed) and the problem in εtatiεtics of finding the "best fit" line to a set of 2-dimensional co-ordinateε which are εubject to experimental error.

The preεent invention uses a cost-function which is, like the sum of squares function, a quadratic function, but differs from the sum of squares function is a significant respect. The present cost- function is a εquare of εum function. Thiε iε defined aε

T J* = J CPe- e + F _e . u) ^a dt 0

The uεe of this second kind of cost-function can lead to simplifications in the solution of the H« optimization problem.

Such cost functions, including in particular the εum of squares function and the square of sum function, may be transformed into the frequency domain. The required minimization iε thereby tranεformed from minimizing a function over a time interval to minimizing a function over frequency. What iε being considered iε the frequency εpectrum or, more preciεely, the εignal εpectral denεitieε which define the power of the signal over frequency.

Such a transformation can conveniently be done by uεing the same general technique aε waε discuεsed earlier for transforming the differential equations describing the behaviour of a system into transfer functions; that is, by uεing the Laplace tranεfer function (if continuouε variableε are being uεed ⁾ or the correεponding discrete operator __~~ (if the εyεtem iε a εampled or diεcrete- interval one). Such a tranεformation converts the above functions into the following forms respectively:

J ₃ = (l/2πj-. j [P _c ..(z- ¹ > + F _c ^(z- ¹ )] dz/z, lzl=l

J _Λ = (l/2πj) j Φ _φtp (z-*> dz/z lzl=l

where Φ (z ^_1 > iε the εpectral denεity of <p(t> = P _c .e(t) + F _c .u(t).

It can be εhown that the minimization of the function J _Λ , which iε the function associated with the present system, can be achieved by solving the two equations in F, G, and H given and discussed above. This iε in contrast to the minimization of the function J ₃ (or J. ), which is a much more difficult problem both in abstract mathematical terms and in practice. *

A cost-function which can be represented by a "sum of squares" function appears at first sight to be seriously deficient. This iε because one way in which the εum P _c (z ^-1 ) e(t> + F _c (z ^_1 > u ⁽ t ⁾ can be minimized iε obviouεly for the two parts P _c (z ^_'l > e(t) and F _e (:_- ^■ > u(t> to be equal and oppoεite, and thiε appears to set no constraint on the valueε of e ⁽ t ⁾ and u(t) individually. Thus the value of the error signal e(t) is apparently uncontrolled, while the whole purpose of the control εystem iε to minimize thiε εignal.

This deficiency is, however, far less significant than would appear from the above consideration, for the following reason. The functions P _c and F _c are, of course, frequency sensitive, and in practice they are chosen so that P _c dominates

at low frequencies and F _e dominates at high frequencies, as indicated informally by Fig. 6. It can thus be seen from the above formula for the power spectrum X that the two terms can interact in the undesirable manner discuεεed above over only a relatively εmall and unimportant part of the whole frequency range.

The general requirement on P _c and F _c iε therefore that P _e should be much larger than F _e at low frequencies and vice versa at high frequencies. A simple example of the manner in which thiε can be achieved iε to take P _c and F _c aε

P _c = <1 - o.z- ¹ >/(l - I.z-'), F _c = β.

Here we can for the moment take α = 0, β = 1, and 1 = 1. (In practice, a func¬ tion for F _c which increases with frequency would be chosen, but the above general condition is satisfied by this very simple form of function for F _c .) The frequency range from 0 to « corresponds to a variation of z from 1 to -.1 anticlockwise around the unit circle from 1 to -1 in the complex plane. A large value of a (i.e. close to 1) results in good εteady state control and εmall control signals but εluggiεh response.

The given form of P _c becomes indefinitely large for zero frequency <z = -1> it can be limited to a large but finite value by taking I aε εlightly less than 1. The parameter I represents integral action; 1 = 1 represents full integral action, 1 = 0 represents none.

In many known H» systems, the control law tends to become constant at high frequencies. The controller may therefore fail to cope adequately with high frequency disturbanceε. The simple choice of F _c above haε thiε disadvantage. However, if the εyεtem being controlled iε such that this disadvantage is sig¬ nificant, then F _c can easily be chosen in the present system to avoid this disadvantage. For example, F _c can be chosen aε

Further, shaping of the response can be achieved. For example, with the ship model diεcuεεed above, it iε desirable for the system to ignore wave action, i.e. for its responεe to disturbances in the frequency band of waves to be reduced. Thiε can be achieved by choosing F _c to be high in that frequency band.

If the εystem is a continuous-time one using the Laplace operator ε, a εtandard form for this is

F _c = (βω _n ^a ε)/(ε ^a + 2ζω _n s + ω„ ² >.

A correεponding form in termε of z~ ^' can be obtained by a conventional transformation from thiε form, e.g. by a bilinear transformation. A similar result can be achieved by choosing P _c to be low for the relevant frequency band.

The preεent εystem in fact resultε in desirable εyεtem properties and a relatively simple control calculation procedure. The system properties are scarcely distinguiεhable from thoεe of a conventional H» control εystem. The control calculation procedure implemented by the preεent εyεtem iε, however, much εi pler than that required by a conventional H» εystem. The present εyεtem iε therefore able to operate adaptively.