WO1999017175A1 - Non-linear dynamic predictive device - Google Patents

Non-linear dynamic predictive device Download PDF

Info

Publication number
WO1999017175A1
WO1999017175A1 PCT/US1998/020295 US9820295W WO9917175A1 WO 1999017175 A1 WO1999017175 A1 WO 1999017175A1 US 9820295 W US9820295 W US 9820295W WO 9917175 A1 WO9917175 A1 WO 9917175A1
Authority
WO
WIPO (PCT)
Prior art keywords
output
input
analyzer
linear
filter
Prior art date
Application number
PCT/US1998/020295
Other languages
French (fr)
Inventor
Hong Zhao
Guillermo Sentoni
John P. Guiver
Original Assignee
Aspen Technology, Inc.
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Aspen Technology, Inc. filed Critical Aspen Technology, Inc.
Publication of WO1999017175A1 publication Critical patent/WO1999017175A1/en

Links

Classifications

    • GPHYSICS
    • G05CONTROLLING; REGULATING
    • G05BCONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/048Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators using a predictor
    • GPHYSICS
    • G05CONTROLLING; REGULATING
    • G05BCONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B17/00Systems involving the use of models or simulators of said systems
    • G05B17/02Systems involving the use of models or simulators of said systems electric

Definitions

  • the present invention pertains to a predictive device that models the dynamic input/output relationships of a physical process, particularly in the process industries such as hydrocarbons, polymers, pulp and paper, and utilities.
  • the predictive device is primarily for multivariable process control, but is also applicable to dynamic process monitoring, or to provide a continuous stream of inferred measurements in place of costly or infrequent laboratory or analyzer measurements.
  • NARX Nonlinear Auto- Regressive with exogenous inputs
  • the process model is identified based on a set of decoupled first order dynamic filters.
  • the use of a group of first order dynamic filters in the input layer of the model enhances noise immunity by eliminating the output interaction found in NARX models.
  • This structure circumvents the difficulty of training a recurrent neural network, while achieving good long-term predictions.
  • using this structure to identify process responses that are second order or higher can result in over sensitive coefficients and in undesirable interactions between the first order filters.
  • this approach usually results in an oversized model structure in order to achieve sufficient accuracy, and the model is not capable of modeling complex dynamics such as oscillatory effects.
  • this first order structure is a special case of a more general nonlinear modeling approach described (Sentoni et al., 1996) that is proven to be able to approximate any discrete, causal, time invari.ant, nonlinear SISO process with fading memory.
  • a Laguerre expansion creates a cascaded configuration of a low pass and several identical band pass first order filters.
  • One of the problems of this approach is that may it require an excessively large degree of expansion to obtain sufficient accuracy. Also, it has not been known until now how to extend this methodology in a practical way to a multi- input system.
  • This invention addresses many essential issues for practical non-linear multivariable MPC. It provides the capability to accurately identify non-linear dynamic processes with a structure that
  • the present invention is a dynamic predictive device that predicts or estimates values of process variables that are dynamically dependent on other measured process variables.
  • This invention is especially suited to application in a model predictive control (MPC) system.
  • the predictive device receives input data under the control of an external device controller.
  • the predictive device operates in either configuration mode or one of three runtime modes -prediction mode, horizon mode, or reverse horizon mode.
  • the primary runtime mode is the prediction mode.
  • the input data are such as might be received from a distributed control system (DCS) as found in a manufacturing process.
  • DCS distributed control system
  • the device controller ensures that a contiguous stream of data from the DCS is provided to the predictive device at a synchronous discrete base sample time.
  • the device controller operates the predictive device once per base sample time and receives the prediction from the output of the predictive device. After the prediction mode output is available, the device controller can switch to horizon mode in the interval before the next base sample time.
  • the predictive device can be operated many times during this interval and thus the device controller can conduct a series of experimental scenarios in which a sequence of input data can be specified by the device controller.
  • the sequence of input data can be thought of as a data path the inputs will follow over a forward horizon.
  • the sequence of predictions at the output of the controller is a predicted output path over a prediction horizon and is passed to the device controller for analysis, optimization, or control.
  • the device controller informs the predictive device at the start of an experimental path and synchronizes the presentation of the path with the operation of the device.
  • horizon mode operates exactly the same way as prediction mode, except that the dynamic states are maintained separately so that the predictive device can resume normal prediction mode operation at the next base sample time.
  • the outputs of the filter units are buffered over the course of the path and .are used during reverse horizon operation of the device.
  • reverse horizon mode The purpose of reverse horizon mode is to obtain the sensitivities of the predictive device to changes in an input path.
  • Reverse horizon mode can only be set after horizon mode operation has occurred.
  • the device controller first informs the predictive device the index of the point in the output path for which sensitivities are required.
  • the device controller then synchronizes the reverse operation of the predictive device with the output of sensitivity data at the input paths of the device.
  • each input is scaled and shaped by a preprocessing unit before being passed to a corresponding delay unit which time-aligns data to resolve dead time effects such as pipeline transport delay.
  • Modeling dead-times is an important issue for an MPC system.
  • prediction horizons are usually set large enough so that both dynamics and dead-time effects are taken into account; otherwise the optimal control path may be based on short term information, and the control behavior may become oscillatory or unstable.
  • the predictive device is predicting a single measurement, .and the dead- time units align data relative to the time of that measurement. If predictions at several measurement points are required, then several predictive devices are used in parallel. During configuration mode, the dead times are automatically estimated using training data collected from the plant.
  • the training method consists of constructing individual auto-regressive models between each input and the output at a variety of dead-times, and choosing the dead time corresponding to the best such model.
  • manual override of the automatic settings is possible and should be used if there is additional process knowledge that allows a more appropriate setting.
  • Each dead time unit feeds a dynamic filter unit.
  • the dynamic filter units are used to represent the dynamic information in the process. Internally the dynamic filter units recursively maintain a vector of states. The states derive their values from states at the previous time step and from the current input value.
  • This general filter type can be represented by what is known to those skilled in the art as a discrete state space equation.
  • the preferred embodiment imposes a much-simplified structure on the filter unit that allows for fast computation for MPC and also allows intelligent override of the automatic settings.
  • This simplified structure is composed of first and second order loosely coupled subfilters, only one of which receives direct input from the corresponding delay unit. The practical identification of this filter structure is an essential part of this invention.
  • the outputs of the dynamic filter units are passed to a non-linear analyzer that embodies a static mapping of the filter states to an output value.
  • the exact nature of the non-linear analyzer is not fundamental to this invention. It can embody a nonlinear mapping such as a Non-linear Partial Least Squares model or a Neural Network, or a hybrid combination of linear model and non-linear model.
  • the preferred embodiment makes use of a hybrid model.
  • the reason for this is that a non- parametric non-linear model identified from dynamic data (such as a neural net) cannot, by its nature, be fully analyzed and validated prior to use.
  • the non-linearity of the model means that different dynamic responses will be seen at different operating points.
  • the process engineer can examine and approve these gains prior to closing the loop on the process and is assured of responses consistent with the true process.
  • the linear dynamic response will be sub-optimal for truly nonlinear processes.
  • the responses of the linear model and the hybrid model can be monitored independently and compared.
  • control can be switched, either automatically or by the operator, back to the safety of the linear model.
  • the output of the non- linear analyzer is passed through a postprocessing unit that converts the internal units to engineering units.
  • a postprocessing unit that converts the internal units to engineering units.
  • This invention is that its structure is shown to be able to approximate a large class of non-linear processes (any discrete, causal, time invariant, nonlinear multi-input/single output (MISO) process with fading memory), but is still simple enough to allow incorporation of process knowledge, is computationally fast enough for practical non-linear MPC, and can be configured with sufficient accuracy in a practical manner.
  • MISO nonlinear multi-input/single output
  • FIG. 7 is an overall block diagram of the invention showing both the runtime and training components.
  • FIG. 2 shows the runtime structure of an individual preprocessing unit.
  • FIG. 3 shows the runtime structure of .an individual delay unit.
  • FIG. 4 shows the forward flow internal decomposition of an individual filter unit into cascaded subfilter units.
  • FIG. 5 shows the preferred forward flow structure of a primary first order subfilter unit.
  • FIG. 6 shows the preferred forward flow structure of a secondary first order subfilter unit and the preferred coupling with the previous subfilter in the cascade.
  • FIG. 7 shows the preferred forward flow structure of a primary second order subfilter unit.
  • FIG. 8 shows the preferred forward flow structure of a secondary second order subfilter unit and the preferred coupling with the previous subfilter in the cascade.
  • FIG. 9 shows a typical feedforward configuration of the non-linear analyzer.
  • FIG.10 shows the reverse flow configuration of the non-linear analyzer depicted in FIG. 9.
  • FIG.11 shows the reverse flow internal decomposition of an individual filter unit into cascaded subfilter units.
  • FIG. 12 shows a method of training an individual delay unit.
  • FIG. 13 shows the first order decoupled structure used at the start of each iteration of the preferred dynamic filter unit identification method.
  • FIG. 14 shows that reverse flow of data through a matrix structure can be described mathematically by forward flow of data through the transpose matrix structure.
  • FIG.l is .an overall block diagram of the invention and its context.
  • An external device controller (50) synchronizes the flow of data to and from the predictive device via the data paths (18), (14), and (64).
  • the device controller also controls the mode of operation and the path stepping of the predictive device via the control path (54).
  • the external device controller may also communicate with a DCS (10) or other data/control system both for requesting data and for requesting control changes to the modeled process; however the exact external context and configuration of the device controller is beyond the scope of this application.
  • index k that represents a data point in a sequence of data points. This index has different meanings depending on whether the forward operational mode of the device is prediction mode or horizon mode.
  • prediction mode data is provided at a regular sampling interval At to the input nodes (18) of the device. Data is passed in a forward direction through the device. For simplicity of notation, the sample point T 0 +kAt is denoted by the index k.
  • horizon mode a sequence of data representing a forward data path is provided to the inputs. This data path may represent a proposed path for manipulated variables for process control purposes, or may represent a holding of the inputs to constant values in order to determine the steady state output of the device. The starting point of this path is taken to be the most recent input sample provided in prediction mode. Index 0 represents this st-arting point and index k represents the k" 1 data point in this path.
  • V.1.1 FORWARD RUNTIME OPERATION OF A PREPROCESSING UNIT Each input feeds a preprocessing unit (20) which is used to convert the engineering units of each data value to a common normalized unit whose lower .and upper limits are, by preference, -1 and 1 respectively, or 0 and 1 respectively.
  • the preprocessing unit can also shape the data by passing it through a nonlinear transformation.
  • u E (k) is the value of an input in engineering units
  • u(k) is the preprocessed value in normalized units.
  • the scale and offset values as stored in the configuration file (30 - FIG. 1) are, in general, different for each input variable, and are determined in the configuration mode.
  • each delayed input value is passed to an individual filter unit (24).
  • the general internal feedforward structure of a filter unit (24) is shown in FIG. 4.
  • the general feedforward structure is composed of S cascaded subfilters (402, 404, ..., 406).
  • the first subfilter in the cascade (400) is referred to as the primary subfilter.
  • Non-primary subfilters are referred to as secondary subfilters. All the subfilters are alike except that the primary subfilter receives no input from another subfilter, and the final subfilter sends no output to another subfilter. Now the general form of the primary subfilter will be described in detail.
  • the primary subfilter maintains a vector (412) of states x t (k) at each time k.
  • An internal single time step delay unit (414) feeds the vector state to a coupling unit (420) and to a matrix unit (416).
  • the matrix unit converts the delayed state vector (418) and feeds it to a vector addition unit (408).
  • the input to the filter unit u d (k) is expanded and linearly scaled by the input coupling unit (410) to a vector of values of the same dimension as the state vector.
  • the vector addition unit then combines its two input streams to produce the vector of states for the current time.
  • the operation just described for the primary subfilter is conveniently described in mathematical matrix and column vector notation as:
  • the subfilter coupling unit (420) determines how state values at time k-1 affect the state units in the next subfilter at time k.
  • the subfilter coupling unit uses the coupling matrix 7 to perform a linear transformation of state vector x,(k-l) which is passed to the vector addition unit of the next subfilter.
  • x s ⁇ k) A s x s (k - l)+ T s x s _,(k - ⁇ ) + b s u d (k) (4)
  • the subfilters are all of first or second order.
  • a first order subfilter maintains just one state.
  • the preferred embodiment for a first order primary subfilter (500) is shown in FIG. 5.
  • the vectorizing unit (502) .and the matrix unit (504) collapse to become scaling operations so that the state vector (506) is represented by:
  • x 1 (k) ⁇ 1 x 1 (k-l)+(l- ⁇ l )u d (k) (5)
  • the preferred embodiment for a first order secondary subfilter (600) is shown in FIG. 6.
  • the secondary subfilter receives no direct input, but instead receives cascaded input from the previous subfilter.
  • the preferred coupling is a loose coupling scheme (602) in which only the last state component of the previous subfilter contributes. Note that the previous subfilter is not required to be a first order subfilter.
  • the state vector (606) is represented by:
  • Second order subfilters maintain two states.
  • the preferred embodiment for a second order primary subfilter (700) is shown in FIG. 7.
  • the state vector x,(k) is shown in terms of its two components x ⁇ (k) (708) andx, 2 (k) (710).
  • the vectorizing unit (702) creates two inputs to the vector addition unit (714), the second of which is fixed at zero.
  • the delayed states (704) and (706) are fed to the matrix unit (712) whose outputs are also fed to the vector addition unit (712) which adds the matrix transformed states to the vectorized inputs producing the current state.
  • the preferred embodiment for a second order secondary subfilter (800) is shown in FIG. 8.
  • the state vector x/k) is shown in terms of its two components x terminat(k) (808) and x l2 (k) (810).
  • the preferred coupling with the previous subfilter unit is a loose coupling scheme (802) in which only the last state component of the previous subfilter contributes to the first state component of the current subfilter. Note that the previous subfilter is not required to be a first order subfilter or second order subfilter.
  • the output of the coupling unit is fed to the addition unit (814).
  • the delayed states (804) and (806) are fed to the state matrix unit (812) whose outputs are also fed to the vector addition unit (812) which adds the state matrix transformed states to the output of the coupling unit, producing the current state. Note that due to the (1,0) structure of the second state matrix row, and the zero second row of the coupling matrix, the current second state component (810) is just equal to the delayed first component (804):
  • the outputs (28) of the filter units (24) provide input to the non-linear analyzer (26).
  • the exact structure and configuration of the non-linear analyzer (26) is not central to this application. It is the interaction of the non-linear analyzer (26) with the filter units (24), and the operation and configuration of the filter units (24) that forms the core of this invention.
  • the preferred embodiment for reasons discussed in the summary of the invention is a hybrid parallel combination of linear and non- linear. However, for clarity of explanation, a standard neural network structure is described which is well known to those skilled in the art. This structure is shown in FIG 9. The equations for this structure are:
  • the postprocessing unit (32) in FIG. 1 is used to scale the output from the normalized units to engineering units.
  • the postprocessing unit can also shape the data by passing it through a non-linear transformation.
  • the preferred embodiment uses a simple scale and offset.
  • the scale and offset represent the mapping from engineering units to normalized units.
  • the reverse horizon mode of operation is only allowed immediately following horizon mode operation.
  • Horizon mode operation buffers the states (28) output by the filter units (24) over the course of the forward path.
  • the purpose of reverse horizon mode is to obtain the sensitivity of any point y(k) of the prediction path (output by the device in horizon mode) with respect to any point in the input path u(l).
  • the mathematical derivatives of the prediction with respect to the inputs are required.
  • the mathematical derivatives measure how sensitive a state is in response to a small change in an input.
  • the dynamic nature of the predictive device means that a change in input at time k will start to have an effect on the output as soon as the minimum dead-time has passed and will continue to have an effect infinitely into the future. In most practical applications systems are identified to have fading memory so that the effect into the future recedes with time.
  • the aim is to plan a sequence of moves for the inputs corresponding to manipulated variables (MVs). The effect of these moves needs to be predicted on the controlled variables (CVs) along a prediction path.
  • MVs manipulated variables
  • a constrained optimization algorithm is then used to find the move sequences that predict an optimal prediction path according to some desired criteria.
  • the external device controller specifies the output path index k.
  • the device then outputs in sequence the sensitivities (64) in reverse order at the input nodes of the device.
  • the sensitivity of the output y E (k) of the device with respect to any variable v is represented by ⁇ . k v .It is this sensitivity value, rather than an external data value that is fed back through the device when operating in reverse horizon mode.
  • the postprocessing unit always receives the value of 1 at its output node in reverse operation.
  • the reverse runtime operation of a neural net model is well known to those skilled in the art and is shown in FIG. 10.
  • the output from the reverse operation of the postprocessing unit ⁇ k y(k) is presented at the output node of the non-linear analyzer (26).
  • ⁇ k ⁇ k ⁇ k ⁇ h ⁇ k)i ⁇ h' ⁇ h ⁇ k))
  • ⁇ k x i (k) ⁇ w hi ⁇ k ⁇ h (k)
  • FIG. 11 shows the filter structure of FIG. 4, but with data flowing in the reverse direction.
  • the vector of sensitivities ⁇ k x s (k) is presented at the output channels (1120, 1122...1124) of the filter unit (24) and cycled in reverse through the filter structure.
  • This reverse operation is indexed by / ⁇ k .
  • the resulting sensitivity ⁇ k u d (l) is output at the input channel (1110) of the filter unit (24).
  • the filter unit (24) itself is not constrained to operate under this assumption.
  • the reverse operation of a delay (1130) is represented by q which is the unit delay in the reverse time direction since the index / is decreasing at each iteration.
  • the reverse operation of a matrix operation (1132, 1134) or a vector operation (1136) is represented mathematically as the transpose of the forward operation.
  • FIG. 14 shows the individual channels represented by a 3x2 -matrix operation which in forward operation maps two input channels to three output channels, and in reverse operation maps three input channels to two output channels.
  • the reverse operation of a delay unit (22) corresponds to a delay in the reverse sequencing:
  • the predictive device is configured, in the preferred embodiment, using training data collected from the process. However, a process engineer can override any automated configuration settings.
  • the training data set should represent one or more data sets which have been collected at the same base-time sample rate that will be used by the external device controller to present data to the predictive device in prediction mode. Each set of data should represent a contiguous sequence of representative.
  • the training of the predictive device is done in stages, each stage representing a major component of the predictive device.
  • the engineering units may be different for each input variable, leading to a different scale and offset for each preprocessing/postprocessing unit.
  • FIG. 12 shows a simple advisory procedure for suggesting delay times. A process engineer can override these advisory settings.
  • d mn and max are user settable limits for the delay time and the procedure calculates a delay time d such that
  • a practical means of configuring a filter unit (24) is an essential aspect of this invention.
  • the preferred method of configuration is initialized using the simplified filter structure shown in FIG. 13 in which all subfilters are first order and decoupled. This is the structure used in (Graettinger, et al, 1994). It is important to note that this structure is used for initialization of the configuration procedure and does not represent the final suggested filter configuration.
  • Step l The operator specifies an appropriate dominant time constant T, associated with each input variable. This can be specified from engineering knowledge or through an automated approach such as Frequency Analysis or a Back Propagation Through Time algorithm.
  • the value of the initial time constant is not critical the proposed configuration method automatically searches the dominant time range for the best values.
  • a fifth order system can be created using:
  • each subfilter (1302, 1304, 1306) yields a corresponding single state (1312, 1314, 1316) which is decoupled from the other subfilter states.
  • This initial filter structure represents the equation
  • Step 4 Balance each subsystem A ( ,b ( ,c ⁇ of the MISO block diagonal system based on controllability & observability theory.
  • the balancing procedure allows order reduction of a state space system by transforming the states so that the controllability and observability properties of the original system are substantially concentrated in the first part of the state vector.
  • Balancing of a linear state space system is a method of reduction well known to those skilled in the art. Other methods of model reduction, such as Hankel reduction, can be substituted. A summary of the balancing method is now given.
  • step 6 If the number of iterations has exceeded the maximum allowable, go to step 6.
  • the eigenvalues 1, • • • ,5 ) of the controllability gramian P, (equivalently the observability gramian) are calculated; these are always positive and real because the controllability gramian is positive definite.
  • the final order S, of each filter unit (24) is then calculated such that
  • is parameter of the training method and is a value less than 1, a good practical value being 0.95.
  • This order represents the total number of states of an individual filter unit (24).
  • the preferred ordering of these subfilter units is according to time-constant, with the fastest unit being the primary subfilter.
  • Another favored approach is to perform model reduction by initializing with Laguerre type filter units as described in section V.4.2, rather than the simple diagonal filter structure of FIG. 13. Sufficient qu.antity of Laguerre type filter units span the full range of dynamics in the process, and thus the iterative process described above is not needed.
  • a non-linear model reduction can be achieved by performing a linear model reduction on the linear system whose states are defined by the Laguerre filters and whose outputs are defined by pre-transformed values at the hidden layer of the neural net: ⁇ ,(k),-, ⁇ H (k)
  • V.3.4 CONFIGURING THE NON-LINEAR ANALYZER The configuration of the non-linear analyzer (26) is not a central aspect of this application.
  • the non-linear analyzer (26) is trained to optimally map the outputs of the filter units (24) to the corresponding t.arget output. Training of a neural net is described in detail in (Bishop, 95) for example.
  • V.4 UNIVERSALITY OF THE PREDICTION DEVICE
  • the predictive device is shown, in this section, to be able to approximate any time invariant, causal, fading memory system (defined below). In order to prove this, some precise notation and definitions will be needed.
  • a variable u represents a vector or a sequence in accordance with the context, while n(k) represents a value of the sequence at the particular time k.
  • R w denotes the normed linear space of real N- vectors (viewed as column vectors) with norm
  • max 1 ⁇ 1 ⁇ JV are denoted in uppercase bold. Functions are denoted in italic lowercase if they are scalars and in bold if they are vector valued.
  • ⁇ (z) (respectively ⁇ (Z + ) .and / ⁇ (Z_)), be the space of bounded Re valued sequences defined on Z (respectively Z + and Z_ ) with the norm:
  • a function F : 1 (Z_ ) -> R is called a functional on / ⁇ (z_ ), and a function 3 : / ⁇ (Z_ ) — > /" (z) is called an operator.
  • a function F : 1 (Z_ ) -> R is called a functional on / ⁇ (z_ )
  • a function 3 : / ⁇ (Z_ ) — > /" (z) is called an operator.
  • Eu rather than E[u] and 3u(&) rather than 3[u](&) .
  • Fading Memory 3 : 1 (z) — » l ⁇ (z) has fading memory on a subset
  • Every sequence u in /" (Z_ )c.an be associated with a causal extension sequence u c in ⁇ (z) defined as: and each time invariant causal operator 3 can be associated with a functional F on / (Z_) defined by
  • the operator 3 can be recovered from its associated functional F via
  • the set of the Laguerre systems is defined in the complex z-transform plane as:
  • L' s (z) is the Z transform of l s ' (k), the s-th order system for the z ' -th input.
  • a, : is the z ' -th input generating pole, such that ⁇ a t I ⁇ 1. This pole is selected as
  • d t is the time delay associated with the z ' -th input variable.
  • ⁇ t : ⁇ -a
  • Each filter unit (24) is configured as a single structured ⁇ A, ,B . ⁇ subfilter.
  • the key point here is that the representation is decoupled by input. Balancing can be done to decrease the order of the Laguerre systems, and similarity transforms can be done on the Laguerre filters in order to simplify the configuration to utilize the preferred subfilter units. Similarity transformations do not affect the accuracy of the representation and so proving that the use of Laguerre filters decoupled by input approximate any time invariant, causal, fading memory system is equivalent to proving the preferred subfilter structure can approximate any such system.
  • the balancing is a practical mechanism to reduce order without degrading performance.
  • E is a compact metric space and G a set of continuous fimctionals on E that separates points, that is for any distinct u, v e E there is a G e G such that
  • Lemma 1 K_ ⁇ ju 6 / (Z_ )
  • K_ is not compact in / ⁇ (Z_) with the usual supremum norm
  • K_ [- /, ⁇ ] the restriction of K__ to [- /, ⁇ ] .
  • K_ [- /, ⁇ ] is uniformly bounded by c x and is composed of a finite set of values, hence compact in / ⁇ [-/, ⁇ ].
  • the G s ' ) separate points in / ⁇ (Z_ ), that is, for any distinct u,v ⁇ (Z_ ) there is a G s ' e G such thatG ⁇ G S ' ⁇ .
  • a time invariant causal operator 3 can be recovered from its associated functional through equation (36) as

Abstract

A non-linear dynamic predictive device (60) is disclosed which operates either in a configuration mode or in one of three runtime modes: prediction mode, horizon mode, or reverse horizon mode. An external device controller (50) sets the mode and determines the data source and the frequency of data. In prediction mode, the input data are such as might be received from a distributed control system (DCS) (10) as found in a manufacturing process; the device controller ensures that a contiguous stream of data from the DCS is provided to the predictive device at a synchronous discrete base sample time. In prediction mode, the device controller operates the predictive device once per base sample time and receives the output from the predictive device through path (14). In horizon mode and reverse horizon mode, the device controller operates the predictive device additionally many times during base sample time interval. In horizon mode, additional data is provided through path (52). In reverse horizon mode data is passed in a reverse direction through the device, utilizing information stored during horizon mode, and returned to the device controller through path (66). In the forward modes, the data are passed to a series of preprocessing units (20) which convert each input variable (18) from engineering units to normalized units.

Description

NON-LINEAR DYNAMIC PREDICTIVE DEVICE
I. FIELD OF THE INVENTION
The present invention pertains to a predictive device that models the dynamic input/output relationships of a physical process, particularly in the process industries such as hydrocarbons, polymers, pulp and paper, and utilities. The predictive device is primarily for multivariable process control, but is also applicable to dynamic process monitoring, or to provide a continuous stream of inferred measurements in place of costly or infrequent laboratory or analyzer measurements.
II. BACKGROUND OF THE INVENTION
Most existing industrial products designed for multivariable model predictive control (MPC) employ linear step-response models or finite impulse response (FIR) models. These approaches result in over-parameterization of the models (Qin and Badgwell, 1996). For example, the dynamics of a first order single input/single output SISO process which can be represented with only three parameters (gain, time constant and dead- time) in a parametric form typically require from 30 to 120 coefficients to describe in a step-response or FIR model. This over-parameterization problem is exacerbated for non-linear models since standard non-parametric approaches, such as Volterra series, lead to an exponential growth in the number of parameters to be identified. An alternative way to overcome these problems for nonlinear systems is the use of parametric models such as input-output Nonlinear Auto- Regressive with exogenous inputs (NARX). Though NARX models are found in many case-studies, a problem with NARX models using feed forward neural networks is that they offer only short-term predictions (Su, et al, 1992). MPC controllers require dynamic models capable of providing long-term predictions. Recurrent neural networks with internal or external feedback connections provide a better solution to the long-term prediction problem, but training such models is very difficult. The approach described in (Graettinger, et al, 1994) and (Zhao, et al, 1997) provides a partial solution to this dilemma. The process model is identified based on a set of decoupled first order dynamic filters. The use of a group of first order dynamic filters in the input layer of the model enhances noise immunity by eliminating the output interaction found in NARX models. This structure circumvents the difficulty of training a recurrent neural network, while achieving good long-term predictions. However, using this structure to identify process responses that are second order or higher can result in over sensitive coefficients and in undesirable interactions between the first order filters. In addition, this approach usually results in an oversized model structure in order to achieve sufficient accuracy, and the model is not capable of modeling complex dynamics such as oscillatory effects. In the single input variable case, this first order structure is a special case of a more general nonlinear modeling approach described (Sentoni et al., 1996) that is proven to be able to approximate any discrete, causal, time invari.ant, nonlinear SISO process with fading memory. In this approach a Laguerre expansion creates a cascaded configuration of a low pass and several identical band pass first order filters. One of the problems of this approach is that may it require an excessively large degree of expansion to obtain sufficient accuracy. Also, it has not been known until now how to extend this methodology in a practical way to a multi- input system.
This invention addresses many essential issues for practical non-linear multivariable MPC. It provides the capability to accurately identify non-linear dynamic processes with a structure that
• has close to minimum parameterization • can be practically identified with sufficient accuracy
• makes good physical sense and allows incorporation of process knowledge
• can be proven to identify a large class of practical processes
• can provide the necessary information for process control III. SUMMARY OF THE INVENTION
The present invention is a dynamic predictive device that predicts or estimates values of process variables that are dynamically dependent on other measured process variables. This invention is especially suited to application in a model predictive control (MPC) system. The predictive device receives input data under the control of an external device controller. The predictive device operates in either configuration mode or one of three runtime modes -prediction mode, horizon mode, or reverse horizon mode.
The primary runtime mode is the prediction mode. In this mode, the input data are such as might be received from a distributed control system (DCS) as found in a manufacturing process. The device controller ensures that a contiguous stream of data from the DCS is provided to the predictive device at a synchronous discrete base sample time. The device controller operates the predictive device once per base sample time and receives the prediction from the output of the predictive device. After the prediction mode output is available, the device controller can switch to horizon mode in the interval before the next base sample time. The predictive device can be operated many times during this interval and thus the device controller can conduct a series of experimental scenarios in which a sequence of input data can be specified by the device controller. The sequence of input data can be thought of as a data path the inputs will follow over a forward horizon. The sequence of predictions at the output of the controller is a predicted output path over a prediction horizon and is passed to the device controller for analysis, optimization, or control. The device controller informs the predictive device at the start of an experimental path and synchronizes the presentation of the path with the operation of the device. Internally, horizon mode operates exactly the same way as prediction mode, except that the dynamic states are maintained separately so that the predictive device can resume normal prediction mode operation at the next base sample time. In addition, the outputs of the filter units are buffered over the course of the path and .are used during reverse horizon operation of the device. The purpose of reverse horizon mode is to obtain the sensitivities of the predictive device to changes in an input path. Reverse horizon mode can only be set after horizon mode operation has occurred. The device controller first informs the predictive device the index of the point in the output path for which sensitivities are required. The device controller then synchronizes the reverse operation of the predictive device with the output of sensitivity data at the input paths of the device.
In forward operation, each input is scaled and shaped by a preprocessing unit before being passed to a corresponding delay unit which time-aligns data to resolve dead time effects such as pipeline transport delay. Modeling dead-times is an important issue for an MPC system. In practical MPC, prediction horizons are usually set large enough so that both dynamics and dead-time effects are taken into account; otherwise the optimal control path may be based on short term information, and the control behavior may become oscillatory or unstable. In the preferred embodiment, the predictive device is predicting a single measurement, .and the dead- time units align data relative to the time of that measurement. If predictions at several measurement points are required, then several predictive devices are used in parallel. During configuration mode, the dead times are automatically estimated using training data collected from the plant. In the preferred embodiment the training method consists of constructing individual auto-regressive models between each input and the output at a variety of dead-times, and choosing the dead time corresponding to the best such model. As with other components of the invention, manual override of the automatic settings is possible and should be used if there is additional process knowledge that allows a more appropriate setting.
Each dead time unit feeds a dynamic filter unit. The dynamic filter units are used to represent the dynamic information in the process. Internally the dynamic filter units recursively maintain a vector of states. The states derive their values from states at the previous time step and from the current input value. This general filter type can be represented by what is known to those skilled in the art as a discrete state space equation. The preferred embodiment imposes a much-simplified structure on the filter unit that allows for fast computation for MPC and also allows intelligent override of the automatic settings. This simplified structure is composed of first and second order loosely coupled subfilters, only one of which receives direct input from the corresponding delay unit. The practical identification of this filter structure is an essential part of this invention. The outputs of the dynamic filter units are passed to a non-linear analyzer that embodies a static mapping of the filter states to an output value. The exact nature of the non-linear analyzer is not fundamental to this invention. It can embody a nonlinear mapping such as a Non-linear Partial Least Squares model or a Neural Network, or a hybrid combination of linear model and non-linear model. The preferred embodiment makes use of a hybrid model. The reason for this is that a non- parametric non-linear model identified from dynamic data (such as a neural net) cannot, by its nature, be fully analyzed and validated prior to use. The non-linearity of the model means that different dynamic responses will be seen at different operating points. If the process being modeled is truly non-linear, these dynamic responses will be an improvement over linear dynamic models in operating regions corresponding to the training data, but may be erroneous in previously unseen operating regions. When the non-linear model is used within the context of MPC, erroneous responses, especially those indicating persistent and invalid gain reversals can create instabilities in the MPC controller. With a hybrid approach, a non-linear model is used to model the errors between the linear dynamic model and the true process. The hybrid dynamic model is a parallel combination of the linear dynamic model with the error correction model. The dynamic response of the linear model can be analyzed completely prior to use, since the gains are fixed and independent of the operating point. The process engineer can examine and approve these gains prior to closing the loop on the process and is assured of responses consistent with the true process. However, the linear dynamic response will be sub-optimal for truly nonlinear processes. In online operation of the hybrid model within an MPC framework, the responses of the linear model and the hybrid model can be monitored independently and compared. In operating regions where the non-linear model shows persistently poor response, control can be switched, either automatically or by the operator, back to the safety of the linear model.
The output of the non- linear analyzer is passed through a postprocessing unit that converts the internal units to engineering units. The importance of this invention is that its structure is shown to be able to approximate a large class of non-linear processes (any discrete, causal, time invariant, nonlinear multi-input/single output (MISO) process with fading memory), but is still simple enough to allow incorporation of process knowledge, is computationally fast enough for practical non-linear MPC, and can be configured with sufficient accuracy in a practical manner.
IV. BRIEF DESCRIPTION OF THE DRAWINGS
The textual description of the present invention makes detailed reference to the following drawings:
FIG. 7 is an overall block diagram of the invention showing both the runtime and training components.
FIG. 2 shows the runtime structure of an individual preprocessing unit.
FIG. 3 shows the runtime structure of .an individual delay unit.
FIG. 4 shows the forward flow internal decomposition of an individual filter unit into cascaded subfilter units. FIG. 5 shows the preferred forward flow structure of a primary first order subfilter unit.
FIG. 6 shows the preferred forward flow structure of a secondary first order subfilter unit and the preferred coupling with the previous subfilter in the cascade.
FIG. 7 shows the preferred forward flow structure of a primary second order subfilter unit.
FIG. 8 shows the preferred forward flow structure of a secondary second order subfilter unit and the preferred coupling with the previous subfilter in the cascade. FIG. 9 shows a typical feedforward configuration of the non-linear analyzer.
FIG.10 shows the reverse flow configuration of the non-linear analyzer depicted in FIG. 9.
FIG.11 shows the reverse flow internal decomposition of an individual filter unit into cascaded subfilter units.
FIG. 12 shows a method of training an individual delay unit.
FIG. 13 shows the first order decoupled structure used at the start of each iteration of the preferred dynamic filter unit identification method.
FIG. 14 shows that reverse flow of data through a matrix structure can be described mathematically by forward flow of data through the transpose matrix structure.
V. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS
FIG.l is .an overall block diagram of the invention and its context. An external device controller (50) synchronizes the flow of data to and from the predictive device via the data paths (18), (14), and (64). The device controller also controls the mode of operation and the path stepping of the predictive device via the control path (54). The external device controller may also communicate with a DCS (10) or other data/control system both for requesting data and for requesting control changes to the modeled process; however the exact external context and configuration of the device controller is beyond the scope of this application.
V.1 FORWARD RUNTIME OPERATION OF THE PREDICTION DEVICE
The figures .and equations in this detailed description refer to an index k that represents a data point in a sequence of data points. This index has different meanings depending on whether the forward operational mode of the device is prediction mode or horizon mode. In prediction mode data is provided at a regular sampling interval At to the input nodes (18) of the device. Data is passed in a forward direction through the device. For simplicity of notation, the sample point T0+kAt is denoted by the index k. In horizon mode, a sequence of data representing a forward data path is provided to the inputs. This data path may represent a proposed path for manipulated variables for process control purposes, or may represent a holding of the inputs to constant values in order to determine the steady state output of the device. The starting point of this path is taken to be the most recent input sample provided in prediction mode. Index 0 represents this st-arting point and index k represents the k"1 data point in this path.
V.1.1 FORWARD RUNTIME OPERATION OF A PREPROCESSING UNIT Each input feeds a preprocessing unit (20) which is used to convert the engineering units of each data value to a common normalized unit whose lower .and upper limits are, by preference, -1 and 1 respectively, or 0 and 1 respectively.
The preprocessing unit can also shape the data by passing it through a nonlinear transformation. However, the preferred embodiment uses a simple scale and offset as shown in FIG. 2 and equation (1): u(k)= suE(k)+ o (1)
where uE(k) is the value of an input in engineering units, and u(k) is the preprocessed value in normalized units. The scale and offset values as stored in the configuration file (30 - FIG. 1) are, in general, different for each input variable, and are determined in the configuration mode.
V.1.2 FORWARD RUNTIME OPERATION OF A DELAY UNIT
Data flows from each preprocessing unit to a corresponding delay unit (22). The forward run-time operation of the delay unit (22) is shown in FIG 3 .and equation (2). The output wrf(^)(304) of an individual delay unit (300) is equal to the input u(k) (302) delayed by d sample times. The value of d may be different for each delay unit (22) and is retrieved from the configuration file (30 - FIG. 1). This may be implemented as a shift register with a tap at the d!h unit. ud(k) = u(k - d) (2) This equation can also be written in terms of the unit delay operator q'\
ud (k) = q~d u(k)
V.1.3 FORWARD RUNTIME OPERATION OF THE FILTER UNITS Referring again to FIG. 1, each delayed input value is passed to an individual filter unit (24). The general internal feedforward structure of a filter unit (24) is shown in FIG. 4. The general feedforward structure is composed of S cascaded subfilters (402, 404, ..., 406). The first subfilter in the cascade (400) is referred to as the primary subfilter. Non-primary subfilters are referred to as secondary subfilters. All the subfilters are alike except that the primary subfilter receives no input from another subfilter, and the final subfilter sends no output to another subfilter. Now the general form of the primary subfilter will be described in detail.
The primary subfilter maintains a vector (412) of states xt(k) at each time k. An internal single time step delay unit (414) feeds the vector state to a coupling unit (420) and to a matrix unit (416). The matrix unit converts the delayed state vector (418) and feeds it to a vector addition unit (408). The input to the filter unit ud(k) is expanded and linearly scaled by the input coupling unit (410) to a vector of values of the same dimension as the state vector. The vector addition unit then combines its two input streams to produce the vector of states for the current time. The operation just described for the primary subfilter is conveniently described in mathematical matrix and column vector notation as:
x1 (k)= Alx1(k -l)+blud{k) (3)
Such an equation is known, to those skilled in the art, as a linear state space equation with a single input. If no structure is imposed on A, or b„ then further subfilters are unnecessary since the cascaded subfilter structure can subsumed into a single complicated primary subfilter. However, the preferred subfilter structures as described below, or similar to those described below, are essential for a practical embodiment and application of the invention. The subfilter coupling unit (420) determines how state values at time k-1 affect the state units in the next subfilter at time k. In mathematical terms, the subfilter coupling unit uses the coupling matrix 7 to perform a linear transformation of state vector x,(k-l) which is passed to the vector addition unit of the next subfilter. The operation of a secondary subfilter is conveniently described in mathematical matrix and vector notations as:
xs {k) = Asxs(k - l)+ Tsxs_,(k - \) + bsud (k) (4)
In the preferred embodiment, the subfilters are all of first or second order. A first order subfilter maintains just one state. The preferred embodiment for a first order primary subfilter (500) is shown in FIG. 5. The vectorizing unit (502) .and the matrix unit (504) collapse to become scaling operations so that the state vector (506) is represented by:
x1(k)= Λ1x1(k-l)+(l-Λl )ud(k) (5)
The preferred embodiment for a first order secondary subfilter (600) is shown in FIG. 6. The secondary subfilter receives no direct input, but instead receives cascaded input from the previous subfilter. The preferred coupling is a loose coupling scheme (602) in which only the last state component of the previous subfilter contributes. Note that the previous subfilter is not required to be a first order subfilter. The state vector (606) is represented by:
Figure imgf000012_0001
where the matrix unit λ s (604) is a scalar.
Second order subfilters maintain two states. The preferred embodiment for a second order primary subfilter (700) is shown in FIG. 7. In this figure, the state vector x,(k) is shown in terms of its two components xπ(k) (708) andx,2(k) (710). The vectorizing unit (702) creates two inputs to the vector addition unit (714), the second of which is fixed at zero. The delayed states (704) and (706) are fed to the matrix unit (712) whose outputs are also fed to the vector addition unit (712) which adds the matrix transformed states to the vectorized inputs producing the current state. Note that due to the (1,0) structure of the second matrix row, and the zero second component of the vectorizing unit component, the current second state component (710) is just equal to the delayed first component (704): xn(k)= anxn (k - \)+ auxl2(k -l)+ (l - an - a12 )ud(k) x k) = xn(k-1)
The preferred embodiment for a second order secondary subfilter (800) is shown in FIG. 8. In this figure, the state vector x/k) is shown in terms of its two components x„(k) (808) and xl2(k) (810). The preferred coupling with the previous subfilter unit is a loose coupling scheme (802) in which only the last state component of the previous subfilter contributes to the first state component of the current subfilter. Note that the previous subfilter is not required to be a first order subfilter or second order subfilter. The output of the coupling unit is fed to the addition unit (814). The delayed states (804) and (806) are fed to the state matrix unit (812) whose outputs are also fed to the vector addition unit (812) which adds the state matrix transformed states to the output of the coupling unit, producing the current state. Note that due to the (1,0) structure of the second state matrix row, and the zero second row of the coupling matrix, the current second state component (810) is just equal to the delayed first component (804):
Figure imgf000013_0001
asiχ si (k - 1)+ s2χ s2 (k ~ l)+ (l " fl.i - as2 k_UflS, (k - 1)
Figure imgf000013_0002
If the device is operating in horizon mode current states along the path are maintained in a separate storage area so as not to corrupt the prediction mode states. In horizon mode, k indexes the input path and the states are initialized at the start of the path (k -0) to the prediction mode states. In addition the states at the output of the filter unit are buffered for use in reverse horizon mode. V.1.4 FORWARD RUNTIME OPERATION OF THE NON-LINEAR ANALYZER Referring again to FIG. 1, the outputs (28) of the filter units (24) provide input to the non-linear analyzer (26). The exact structure and configuration of the non-linear analyzer (26) is not central to this application. It is the interaction of the non-linear analyzer (26) with the filter units (24), and the operation and configuration of the filter units (24) that forms the core of this invention. The preferred embodiment, for reasons discussed in the summary of the invention is a hybrid parallel combination of linear and non- linear. However, for clarity of explanation, a standard neural network structure is described which is well known to those skilled in the art. This structure is shown in FIG 9. The equations for this structure are:
JV ξh{k) = who + ∑ Whixi{k) i=l ηh{k)= Xxύι{ξh{k)) (9)
y{k)= ∑J c^k{k) h=\
V.1.5 FORWARD RUNTIME OPERATION OF THE POSTPROCESSING UNIT The postprocessing unit (32) in FIG. 1 is used to scale the output from the normalized units to engineering units. The postprocessing unit can also shape the data by passing it through a non-linear transformation. However, the preferred embodiment uses a simple scale and offset. For consistency with the preprocessing units, the scale and offset represent the mapping from engineering units to normalized units.
yE(k)= Ly(k)-- (10)
The scale and offset values as stored in the configuration file (30 - FIG. 1) and are determined in the configuration mode. V.2 REVERSE RUNTIME OPERATION OF THE PREDICTION DEVICE
The reverse horizon mode of operation is only allowed immediately following horizon mode operation. Horizon mode operation buffers the states (28) output by the filter units (24) over the course of the forward path. The purpose of reverse horizon mode is to obtain the sensitivity of any point y(k) of the prediction path (output by the device in horizon mode) with respect to any point in the input path u(l).
In order to use the invention for process control applications, the mathematical derivatives of the prediction with respect to the inputs are required. The mathematical derivatives measure how sensitive a state is in response to a small change in an input. The dynamic nature of the predictive device means that a change in input at time k will start to have an effect on the output as soon as the minimum dead-time has passed and will continue to have an effect infinitely into the future. In most practical applications systems are identified to have fading memory so that the effect into the future recedes with time. For MPC applications the aim is to plan a sequence of moves for the inputs corresponding to manipulated variables (MVs). The effect of these moves needs to be predicted on the controlled variables (CVs) along a prediction path. A constrained optimization algorithm is then used to find the move sequences that predict an optimal prediction path according to some desired criteria. In reverse horizon mode, the external device controller specifies the output path index k. The device then outputs in sequence the sensitivities (64) in reverse order at the input nodes of the device. In the detailed description below, the sensitivity of the output yE(k) of the device with respect to any variable v is represented by Ω.kv .It is this sensitivity value, rather than an external data value that is fed back through the device when operating in reverse horizon mode.
V.2.1 REVERSE RUNTIME OPERATION OF THE POSTPROCESSING UNIT
The reverse operation of the postprocessing unit (32) is to scale data received at its output node using the inverse of the feedforward scaling shown in equation (10): Ωky(k) = sΩkyE(k) (11)
Since the sensitivity of the output with respect to itself is: ΩkyE{k)= l (12)
the postprocessing unit always receives the value of 1 at its output node in reverse operation.
V.2.2 REVERSE RUNTIME OPERATION OF THE NON-LINEAR ANALYZER
The reverse runtime operation of a neural net model is well known to those skilled in the art and is shown in FIG. 10. The output from the reverse operation of the postprocessing unit Ωk y(k) is presented at the output node of the non-linear analyzer (26). The information flows in a reverse manner through the non-linear analyzer (26) and the resulting sensitivities (62) are output at the input nodes of the non-linear analyzer (26):
Figure imgf000016_0001
Ωkξ k) = Ωkηh{k)i∞h'{ξh{k))
Figure imgf000016_0002
Ωkxi(k)= ∑whiΩkξh(k)
A=l
V.2.3 REVERSE RUNTIME OPERATION OF A FILTER UNIT The effect of a change in the delayed input ud(l) on a the sequence of states being output from a filter unit (24) in horizon mode is complex due to the dependencies of a subfilter's states based on the previous subfilter's states and on the subfilter's previous states. An efficient solution can be derived using the chain rule for ordered derivatives (Werbos, 1994) and is achieved by the reverse operation of the filter unit (24). In reverse horizon mode, the output of each filter unit (24) receives the vector of sensitivities Ωsxs (k) propagated back from the non-linear analyzer (26) operating in reverse mode: (14)
Figure imgf000017_0001
The operation of these equations is shown in FIG. 11, which shows the filter structure of FIG. 4, but with data flowing in the reverse direction. Given the point k in the output path for which the sensitivities are being calculated, the vector of sensitivities Ωkxs (k) is presented at the output channels (1120, 1122...1124) of the filter unit (24) and cycled in reverse through the filter structure. This reverse operation is indexed by / < k . At each iteration / , the resulting sensitivity Ωkud(l) is output at the input channel (1110) of the filter unit (24). For / < k the external input at the output channels (1120, 1122...1124) is in practice zero vector since Ωkxs (/) = 0 . However, the filter unit (24) itself is not constrained to operate under this assumption.
In FIG. 11, the reverse operation of a delay (1130) is represented by q which is the unit delay in the reverse time direction since the index / is decreasing at each iteration. The reverse operation of a matrix operation (1132, 1134) or a vector operation (1136) is represented mathematically as the transpose of the forward operation. The physical justification for this is shown in FIG. 14 which shows the individual channels represented by a 3x2 -matrix operation which in forward operation maps two input channels to three output channels, and in reverse operation maps three input channels to two output channels. V.2.4 REVERSE RUNTIME OPERATION OF A DELAY UNIT
The reverse operation of a delay unit (22) corresponds to a delay in the reverse sequencing:
Ωku(l) = Ωkud (l + d) (15)
V.2.5 REVERSE RUNTIME OPERATION OF A PREPROCESSING UNIT
The reverse operation of a preprocessing unit (20) is to scale data received at its output node using the inverse of the feedforward scaling shown in equation (1):
ΩkuE(!) = -Ωku(l) (16) s
V.3 CONFIGURATION MODE The predictive device is configured, in the preferred embodiment, using training data collected from the process. However, a process engineer can override any automated configuration settings. The training data set should represent one or more data sets which have been collected at the same base-time sample rate that will be used by the external device controller to present data to the predictive device in prediction mode. Each set of data should represent a contiguous sequence of representative.
In order to allow operator approval or override of the configuration settings, the training of the predictive device is done in stages, each stage representing a major component of the predictive device.
V.3.1 CONFIGURING THE PREPROCESSING AND POSTPROCESSING UNITS The scale and offset of a preprocessing or postprocessing unit is determined from the desire to map the minimum Emin and maximum Emax of the corresponding variable's engineering units to the minimum Nmin and maximum Nmax of the normalized units: „ N^ max - N nun
E max - F nun (Λ H\
_ E maxN„ rm.n -E mm N„ max
F max - F nun
The preferred normalized units have Nmm=-\, Nmax=\ . The engineering units may be different for each input variable, leading to a different scale and offset for each preprocessing/postprocessing unit.
V.3.2 CONFIGURING A DELAY UNIT
The configuration of a delay unit (22) is not a central aspect of this application. FIG. 12 shows a simple advisory procedure for suggesting delay times. A process engineer can override these advisory settings. In this procedure dmn and max are user settable limits for the delay time and the procedure calculates a delay time d such that
V.3.3 CONFIGURING A FILTER UNIT
A practical means of configuring a filter unit (24) is an essential aspect of this invention. The preferred method of configuration is initialized using the simplified filter structure shown in FIG. 13 in which all subfilters are first order and decoupled. This is the structure used in (Graettinger, et al, 1994). It is important to note that this structure is used for initialization of the configuration procedure and does not represent the final suggested filter configuration.
Step l The operator specifies an appropriate dominant time constant T, associated with each input variable. This can be specified from engineering knowledge or through an automated approach such as Frequency Analysis or a Back Propagation Through Time algorithm. The value of the initial time constant is not critical the proposed configuration method automatically searches the dominant time range for the best values. Ste 2
For each input, initialize the filter structure in FIG. 13 using a high order system where a number of first order filters are created around the given dominant time constant (dominant frequency, dominant dynamics). For example, a fifth order system can be created using:
Ar λiX= e 0.57;
At λi2= e 0.757;
Δ? λn= e T, (18)
At λu= e 1.257;
M 1 „ L57.
In this simple filter structure, each subfilter (1302, 1304, 1306) yields a corresponding single state (1312, 1314, 1316) which is decoupled from the other subfilter states. This initial filter structure represents the equation
x(k)= Ax(k -l) + ud (k) (19)
which has a simplified diagonal block structure of the form
Figure imgf000020_0001
where
Figure imgf000021_0001
λ. 0 0 0
0 λΛ 0
(21)
0 0
0 0 λ.
1 -JL,
\ -λ, b, =
\ -λ.
Step 3
Map the contiguous input training data through the delay units (22) and filter structure (24) to obtain a set of training state vectors { (&)|& = 1, • • • ,T . Then find a vector c that provides the best linear mapping of the states to the corresponding target outputs ι (&)|& = 1, • • • , T) . One way of doing this is to use the Partial Least Squares method that is well known to those skilled in the art. This results in a multi- input, single-output (MISO) state space system |A,b,cr | in which equations (19), (20), and (21) are supplemented by the equation:
y(k) = cτx (22)
where
Figure imgf000021_0002
Step 4 Balance each subsystem A( ,b( ,c }of the MISO block diagonal system based on controllability & observability theory. The balancing procedure allows order reduction of a state space system by transforming the states so that the controllability and observability properties of the original system are substantially concentrated in the first part of the state vector.
For each input variable, indexed by i, perform the balancing procedure on the sub-system A(,br,c }. Balancing of a linear state space system is a method of reduction well known to those skilled in the art. Other methods of model reduction, such as Hankel reduction, can be substituted. A summary of the balancing method is now given.
For each sub-system JA, ,b( ,c }, compute the controllability and observability Gramians P, > 0, Q, > 0 that satisfy the equations:
Figure imgf000022_0001
Find a matrix R,, using the Cholesky factorization method, such that
Figure imgf000022_0002
Using the singular value decomposition method, diagonalize to obtain the following decomposition:
Figure imgf000022_0003
Define
Figure imgf000022_0004
then
T,P,TX (τ,r)"1Q,TX ∑, (28)
.and the balanced subsystem is obtained through a similarity transform on the states as:
Figure imgf000023_0001
Step 5
Using balanced subsystems find out dominant time constant for each input by reducing each balanced model to a first order model. This is done by considering the dynamics of all but the first state of each input's filter unit (24) to have reached steady state. This leads to:
At
T. = (30)
Ha,) where
a, = α,„ + a ^ll ^ - ll ) */21 (31)
and
Figure imgf000023_0002
Check the convergence of the dominant time constant estimation: If
current previous \ ^
(33)
or the number of iterations has exceeded the maximum allowable, go to step 6.
Otherwise, return to step 2. The maximum number of iterations and ε are parameters of the training method.
Step 6
Once an accurate estimate of the dominant time constant is available for each input variable, the eigenvalues
Figure imgf000023_0003
= 1, • • • ,5 ) of the controllability gramian P, (equivalently the observability gramian) are calculated; these are always positive and real because the controllability gramian is positive definite. The final order S, of each filter unit (24) is then calculated such that
Figure imgf000024_0001
where θ is parameter of the training method and is a value less than 1, a good practical value being 0.95. This order represents the total number of states of an individual filter unit (24).
After determining the model order, truncate the A, matrix so that just the first S, states are used; this truncation is done by selecting the upper left S,x S, submatrix of A, . Then calculate the S, eigenvalues of the truncated A, matrix
Figure imgf000024_0002
= 1, • • • , S, } . Now configure each filter unit (24) using the preferred first and second order subfilter configurations with the preferred couplings as shown in FIG 5 through FIG. 8. Use a first order filter for each real eigenvalue. Use a second order filter for each pair of complex eigenvalues
Figure imgf000024_0003
, where, in FIG. 7 (equation 7) or FIG. 8 (equation 8):
a,. = Λ+ λ " _ (35) al2 = -λλ
The preferred ordering of these subfilter units is according to time-constant, with the fastest unit being the primary subfilter.
Another favored approach is to perform model reduction by initializing with Laguerre type filter units as described in section V.4.2, rather than the simple diagonal filter structure of FIG. 13. Sufficient qu.antity of Laguerre type filter units span the full range of dynamics in the process, and thus the iterative process described above is not needed. In fact a non-linear model reduction can be achieved by performing a linear model reduction on the linear system whose states are defined by the Laguerre filters and whose outputs are defined by pre-transformed values at the hidden layer of the neural net: ξ ,(k),-,ξ H(k)
V.3.4 CONFIGURING THE NON-LINEAR ANALYZER The configuration of the non-linear analyzer (26) is not a central aspect of this application. The non-linear analyzer (26) is trained to optimally map the outputs of the filter units (24) to the corresponding t.arget output. Training of a neural net is described in detail in (Bishop, 95) for example.
V.4 UNIVERSALITY OF THE PREDICTION DEVICE The predictive device is shown, in this section, to be able to approximate any time invariant, causal, fading memory system (defined below). In order to prove this, some precise notation and definitions will be needed.
V.4.1 NOTATION AND DEFINITIONS FOR UNIVERSALITY PROOF
Let Z denote the integers, Z+ the non-negative integers and Z_ the non- positive integers respectively. A variable u represents a vector or a sequence in accordance with the context, while n(k) represents a value of the sequence at the particular time k.
For any positive integer ? > 0, Rw denotes the normed linear space of real N- vectors (viewed as column vectors) with norm |u| = max1≤1≤JV
Figure imgf000025_0001
are denoted in uppercase bold. Functions are denoted in italic lowercase if they are scalars and in bold if they are vector valued.
Let ^(z) (respectively ^(Z+) .and /^(Z_)), be the space of bounded Re valued sequences defined on Z (respectively Z+ and Z_ ) with the norm:
Figure imgf000025_0002
For every decreasing sequence w : Z+ -» (θ,l], lim w(&) = 0 define the following weighted norm:
IH = suPtez.lιι(ArM-*)
A function F : 1 (Z_ ) -> R is called a functional on /^ (z_ ), and a function 3 : /^ (Z_ ) — > /" (z) is called an operator. As a notational simplification the parentheses around the arguments of fimctionals and operators are usually dropped; for example, Eu rather than E[u] and 3u(&) rather than 3[u](&) .
Two specific operators are important. The delay operator defined by
Qdu(k)≡ u(k - d)
and the truncation operator defined by
, , fu(ifc) k ≤ O Pu(k) ≡ { ' [0 k > 0
The following definitions make precise the terms used to characterize the class of systems approximated by the predictive device.
Time invariant: An operator 3 is time-invariant if Qd3 = 3Qd W e Z .
Causality: 3 is causal if u(/) = v(/)V/ < k => 3u(&) = 3v(&) .
Fading Memory: 3 : 1 (z) — » l (z) has fading memory on a subset
K_ /^ (Z_ ) if there is a decreasing sequence w : Z+ — (θ,l], lim w(k) = 0 , such that k→∞ for each u, v e K_ and given ε > 0 there is a δ > 0 such that
|u(jfc)- v(ikjw < £ => |3 (θ)- 3v(θ)| < <?
Every sequence u in /" (Z_ )c.an be associated with a causal extension sequence uc in ^(z) defined as:
Figure imgf000027_0001
and each time invariant causal operator 3 can be associated with a functional F on / (Z_) defined by
Eu = 3uc(θ)
The operator 3 can be recovered from its associated functional F via
3u(*) = £r*u (36)
Then, 3 is continuous if and only if F is, so the above equations establish a one to one correspondence between time invariant causal continuous operators and fimctionals F on /^(z_). In the next the definition of the Laguerre system is given. These can be configured in the general filter structure of FIG. 4 but also have important theoretical properties.
V.4.2 LAGUERRE SYSTEMS
The set of the Laguerre systems is defined in the complex z-transform plane as:
η, z -</,+!
X ~ aj Z
L = s = 0X — ,∞, i = l, — ,N z - a .
where:
L's (z) : is the Z transform of ls' (k), the s-th order system for the z'-th input.
a, : is the z'-th input generating pole, such that \at I < 1. This pole is selected as
at = 1 , where Tf is the dominant time constant for the z'-th input
variable. dt : is the time delay associated with the z'-th input variable. ηt : = \-a)
The whole set of Laguerre systems can be expressed in a state space form that shows a decoupled input form and therefore can be mapped to the general filter structure in FIG 4. Each filter unit (24) is configured as a single structured {A, ,B .} subfilter. The structure of A,, is a lower triangular matrix, and b, = [l 0 • • θ]τ .
The key point here is that the representation is decoupled by input. Balancing can be done to decrease the order of the Laguerre systems, and similarity transforms can be done on the Laguerre filters in order to simplify the configuration to utilize the preferred subfilter units. Similarity transformations do not affect the accuracy of the representation and so proving that the use of Laguerre filters decoupled by input approximate any time invariant, causal, fading memory system is equivalent to proving the preferred subfilter structure can approximate any such system. The balancing is a practical mechanism to reduce order without degrading performance.
V.4.3 PROOF OF APPROXIMATION ABILITY OF LAGUERRE SYSTEMS First some preliminary results are stated:
Stone- eierstrass Theorem (Bovd.1985).
Suppose E is a compact metric space and G a set of continuous fimctionals on E that separates points, that is for any distinct u, v e E there is a G e G such that
Gu ≠ Gx . Then for any continuous functional F on E and given ε > 0 , there are
( i N fimctionals, {G,1 , • • • Gs l > , • • • , G" , ■ ■ ■ G"N } c G , S = ∑ S, and a polynomial
;=1 p : Rs - R , such that for all u € E
Figure imgf000028_0001
The reason for the group indexing, which is not necessary for a general statement of the Stone- Weierstrass theorem, will become apparent in Lemma 2 when each block with a Laguerre operator. In addition, three lemmas are necessary before the theorem can be proved.
Lemma 1 : K_ ≡ ju 6 / (Z_ )|θ < |u| < c, ), is compact with the ||w norm.
Proof: Let u^ be any sequence in K_ . We will find a u(0) e K_ and a subsequence
of p' converging in the \\\\w norm to u °\ It is well know that K_ is not compact in /^(Z_) with the usual supremum norm || (Kolmogorov, 1980). For each /, let be K_ [- /,θ] the restriction of K__ to [- /,θ] . K_ [- /,θ] is uniformly bounded by cx and is composed of a finite set of values, hence compact in /^[-/,θ]. Since ^_[-/,θ] is compact for every /, we can find a subsequence u^ of v p' and a u^0' e 7T_[-/,θ] along which u^ converges:
sup u p (k)- u{o)(k] → 0 as m → ∞ (37)
-l≤k≤O
Now, let ε > 0. Since w(k) -» 0 as k -> oo , we cian find m0 > 0 a such that w(m0 ) ≤ ε/cλ . Since u^"',u^ € K_ , we have that
sup |u(,, (&)- u{0)(k w(- k) ≤ 2c, w(m0 ) < ε (38) k≤-m0
Now from equation (37) we can find mx such that
sup U (p (*)-u(0)(ifc)| < * for m > m (39)
-m0 <jt≤0:
so by equation (38) and equation (39) we can conclude that
Figure imgf000029_0001
which proves that K_ is compact. Lemma 2. The set of functional Gs' ) associated to the discrete Laguerre Operators are continuous with respect to || norm, that is, given any ε > 0 there exists δ > 0 such that llu - v < δ G' -Glv \<ε
Proof. Consider the functional Gs' (•) associated with the Laguerre operator s[-j . Given ε > 0 , chose δ > 0 such that:
\u. -v, <δ G'u -G' \<ε (40)
This is possible due to the continuity of the one dimensional Laguerre operators with respect to the weighted norm as shown in (Sentoni et al, 1996). Therefore, from equation (40) and the definition of the fimctionals
u-v <δ=>\u -v,\ <δ=> G;U-G;V| = |GX-G;V,|< (4i)
which proves Lemma 2
Lemma 3. The Gs' ) separate points in /^ (Z_ ), that is, for any distinct u,vε (Z_ ) there is a Gs' e G such thatG ≠GS'\ .
Proof. Suppose u,ve (Z_ ) are equal except for the z'-th component. Then
Figure imgf000030_0001
by the definition of the fimctionals. It is known from one dimensional theory (Sentoni et al, 1996) that for any distinct w,,v, e /M(Z_) there is a Gs' such that
Gs'ul ≠ Gs'vl ; this result together with equation (42) proves Lemma 3
Approximation Theorem Now given ε > 0, Lemmas 1, 2, 3 together with the Stone- Weierstrass theorem imply that given any continuous functional F on K_ , there is a polynomial p:Rs →R, such that for all UGK_
Fu-p(Gl 1u,--,Glu,-,Gl Nu,--,G ]<ε (43) Because the Laguerre systems are continuous and acting on a bounded space, the G^u are bounded real intervals on so the polynomial * can be replaced by any static model that acts as a universal approximator on a bounded input space, for example, a neural net. In other words (43) can be replaced by
|Eu - NN(G1 ,u,- -, G iU,- - -,G1V- - -,Gs ] < <5- (44)
A time invariant causal operator 3 can be recovered from its associated functional through equation (36) as
3u(k) = FPQ-ku
Now let u € K and k e Z , so PQ'ku e K_ , hence
Figure imgf000031_0001
Since the last equation is true for all k e Z , we conclude that for all u e K_
3u - 3u < ε
In other words, it is possible to approximate any nonlinear discrete time invariant operator having fading memory on K , with a finite set of discrete Laguerre systems followed by a single hidden layer neural net. This completes the proof.
V.5 EQUIVALENTS
Although the foregoing details refer to particular preferred embodiments of the invention, it should be understood that the invention is not limited to these details. Substitutions and alterations, which will occur to those of ordin.ary skill in the art, can be made to the detailed embodiments without departing from the spirit of the invention. These modifications are intended to be within the scope of the present invention.

Claims

CLAIMS What is claimed is:
1. A predictive device for modeling a non-linear, causal, multiple-input single- output system or process, comprising: a plurality of preprocessing units for receiving a working signal including control data inputs, the preprocessing units normalizing the control data inputs, resulting in preprocessed inputs; a plurality of delay units coupled to the preprocessing units, the delay units time aligning the preprocessed inputs, resulting in time aligned inputs; a plurality of filter units coupled to the delay units, the filter units being of a substantially simplified configuration as compared to a configuration based upon discrete state space equations, the filter units filtering the time aligned inputs at least according to time, resulting in filtered states; a non-linear analyzer coupled to the filter units and accepting the filtered states, the non-linear analyzer generating a single analyzer output; a postprocessing unit coupled to the non-linear analyzer to receive the generated analyzer output, the postprocessing unit converting the single analyzer output to a single device output that represents an estimate or prediction of the output of the multiple-input single-output dynamic system being modeled by the device, and wherein the predictive device operates in a plurality of selectable modes including a configuration mode and multiple runtime modes.
2. The device of Claim 1 wherein: the data generated by the predictive device, in any of the selectable modes, are received by a device controller for analysis, monitoring, optimization or control of the modeled process.
3. The device of Claim 1 wherein the preprocessing units normalize the control data inputs by scaling and offsetting the control data inputs, resulting in preprocessed inputs.
4. The device of Claim 1 wherein the postprocessing unit normalizes the analyzer output by scaling and offsetting the analyzer output, resulting in a postprocessed device output .
5. The device of Claim 1 further comprising a plurality of training units for configuring the predictive device in configuration mode, the training units including: a preprocessing training unit to set overrideable parameters, including scale and offset settings, in the plurality of preprocessing units; a delay training unit to set overrideable delay times to the plurality of delay units; a filter training unit to configure the plurality of filter units; a non-linear analyzer training unit to train the non-linear analyzer to optimally map the filtered states to the analyzer output; a postprocessing training unit to set overrideable parameters, including scale and offset settings, in the plurality of postprocessing units; and wherein the training units are activated when the predictive device is operated in configuration mode.
6. The device of Claim 5 wherein the filter training unit is used to configure, in a practical manner allowing relatively simple override by an operator, the plurality of filter units by using input training data for control data input, and by using a corresponding measured process output for target device output training data, and by:
(a) setting an initial estimated dominant time constant associated with each input training data;
(b) for each input training data, initializing a filter structure based on the input training data's initial estimated dominant time constant; (c) for each target device output training data, determining a corresponding target analyzer output training data
(d) mapping input training data through the delay units and the filter structure, resulting in a plurality of training filter state vectors, creating a vector trained to linearly map training states to the corresponding target analyzer outputs, resulting in a multi -input, single-output state space system consisting of a block structure of single-input single-output subsystems, one for each filter unit;
(e) applying a method of order reduction to each single-input, single-output sub-system, resulting in a set of reduced order single-input, single-output sub-systems;
(f) reducing each reduced single-input, single-output sub-system to a first order system to determine an updated estimated dominant time constant for each input training data;
(g) repeating steps (a) through (f) for one or more iterations using the updated estimated dominant time constants in place of the initial estimated dominant time constants; and
(h) configuring each filter unit using the corresponding reduced order single- input, single-output sub-system.
7. The device of Claim 5 wherein the filter training unit is used to configure, in a practical manner allowing relatively simple override by an operator, the plurality of filter units by using input training data for control data input, and by using a corresponding measured process output for target device output training data, and by:
(a) setting an estimated dominant time constant associated with each input training data;
(b) for each input training data, initializing a filter structure based upon the input training data's estimated dominant time constant; (c) for each target device output training data, determining a corresponding target analyzer output training data
(d) mapping input training data through the delay units and the filter structure, resulting in a plurality of training filter state vectors, creating a linear or non-linear analyzer trained to map training states to the corresponding target analyzer outputs, where said analyzer has a set of internal states which are linearly dependent on a current filter state vector, resulting in a linear multi-input, multi-output state space system mapping input training data to said internal analyzer states, where said linear multi- input, multi-output state space system consists of a block structure of single-input multi-output sub-systems, one for each filter unit;
(e) applying a method of order reduction to each single-input, multi-output state space system, resulting in a set of reduced order single-input, multi- output sub-systems; and
(f) configuring each filter unit using the corresponding reduced order single- input, multi-output sub-systems.
8. The device of claim 7 wherein the non- linear analyzer is a neural network and the internal states are summation values of the neural network hidden layer.
9. The device of claim 7 wherein the linear or non-linear analyzer is a Partial Least Squares model and the internal states are latent variable scores.
10. The device of Claim 6 or Claim 7 wherein the method of reduction is a Hankel model reduction procedure.
11. The device of Claim 6 or Claim 7 wherein the method of reduction is an Internally Balanced model reduction procedure.
12. The device of Claim 6 or Claim 7 wherein the method of reduction is an iterative application of any of a variety of model reduction procedures including a Hankel model reduction procedure and an Internally Balanced model reduction procedure.
13. The device of Claim 1 wherein the plurality of selectable runtime modes includes a predictive mode in which:
(i) the predictive device receives a contiguous stream of control data inputs at asynchronous discrete base sample time; and
(ii) the predictive device is operated once per base sample time.
14. The device of claim 13 wherein the contiguous stream of control data inputs is passed from a device controller and the analyzer output is received by the device controller for analysis, monitoring, optimization or control of the modeled process.
15. The device of Claim 1 wherein the plurality of selectable runtime modes comprises an horizon mode in which the predictive device: receives an externally defined sequence of trial future data inputs proceeding from a current prediction mode device state; is operated in response to this trial sequence of data inputs producing a corresponding sequence of at least filtered states, and possible other state information; and stores the filtered states and other state information for use in reverse horizon mode.
16. The device of claim 15 wherein the horizon mode is run one or more times between runs of the predictive device in the predictive mode.
17. The device of claim 15 wherein a contiguous stream of external trial data inputs is passed to the predictive device from a device controller; and the predictions generated during horizon mode are received by the device controller for analysis, monitoring, optimization or control of the modeled process.
18. The device of Claim 1 wherein the plurality of selectable runtime modes comprises a reverse horizon mode in which the predictive device uses (i) the filtered states and other state information from the most recent horizon mode run, and
(ii) an output path index indicating a point in a generated sequence of predictions to obtain the sensitivities of the predictive device to changes in the trial input data sequence used by a most recent horizon mode run, based upon running the predictive device backwards.
19. The device of claim 18 wherein the reverse horizon mode is run one or more times between runs of the predictive device in the predictive mode.
20. The device of claim 18 wherein the predictive device sensitivities generated during reverse horizon mode are received by a device controller for analysis, monitoring, optimization or control of the modeled process.
21. The device of claim 18 wherein a device controller specifies the output path index.
22. The device of Claim 6 or Claim 7 wherein the initializing of the filter structure uses a Laguerre expansion.
23. The device of Claim 1 wherein the plurality of filter units comprise: first and/or second order subfilters.
24. The device of Claim 1 wherein the non-linear analyzer comprises: a neural network.
25. The device of Claim 1 wherein the non-linear analyzer comprises: a linear or non-linear Partial Least Squares model.
26. The device of Claim 1 wherein the non-linear analyzer comprises: a hybrid parallel combination of a linear model with a non-linear model.
27. A computer method for modeling a non-linear, causal, multiple-input single- output, system or process, comprising the steps of: (a) receiving and normalizing a working signal including control data inputs, resulting in preprocessed inputs;
(b) aligning the preprocessed inputs, resulting in time aligned inputs;
(c) using a plurality of filter units, filtering the time aligned inputs, at least according to time, resulting in filtered states;
(d) generating an analyzer output based upon the filtered states, said generating employing a non-linear analyzer; and
(e) converting the analyzer output to a model output that represents an estimate or prediction of the output of the multiple-input single-output dynamic system being modeled by the method.
28. The method of Claim 27 wherein the step of receiving includes receiving a contiguous stream of control data inputs from an external system, said data inputs representing measurements from the modeled process; and further comprising the step of passing model output to an external device or method for analysis, monitoring, optimization or control of the modeled process.
29. The method of Claim 27 wherein the normalizing step employs a scale and offset for each input.
30. The method of Claim 27 wherein the converting step employs a scale and offset.
31. The method of Claim 27 further comprising the step of configuring including the steps of:
(a) setting overrideable parameters, including scale and offset settings, for the receiving and normalizing step; (b) setting overrideable delay times, for the aligning the preprocessed inputs step;
(c) training the plurality of filter units used in the filtering of the time aligned inputs step; (d) training the non-linear analyzer to optimally map the filtered states to the analyzer output in the generating analyzer outputs step; and
(e) setting overrideable parameters, including scale and offset settings, for the converting the analyzer output to a model output.
32. The method of Claim 31 wherein the step of training the filter units includes: employing a filter training unit to configure, in a practical manner allowing relatively simple override by an operator, the plurality of filter units by using input training data for control data input, and by using a corresponding measured process output for target model output training data, and by:
(a) setting an initial estimated dominant time constant associated with each input training data;
(b) for each input training data, initializing a filter structure based on the input training data's initial estimated dominant time constant; (c) for each target device output training data, determining a corresponding target analyzer output training data
(d) mapping input training data through the delay units and the filter structure, resulting in a plurality of training filter state vectors, creating a vector trained to linearly map training states to the corresponding target analyzer outputs, resulting in a multi-input, single-output state space system consisting of a block structure of single-input single-output subsystems, one for each filter unit; (e) applying a method of order reduction to each single-input, single-output sub-system, resulting in a set of reduced order single-input, single-output sub-systems;
(f) reducing each reduced single-input, single-output sub-system to a first order system to determine an updated estimated dominant time constant for each input training data;
(g) repeating steps (a) through (f) for one or more iterations using the updated estimated dominant time constants in place of the initial estimated dominant time constants; and (h) configuring each filter unit using the corresponding reduced order single- input, single-output sub-system.
33. The method of Claim 31 wherein the step of training the filter units includes: employing a filter training unit to configure, in a practical manner allowing relatively simple override by an operator, the plurality of filter units by using input training data for control data input, and by using a corresponding measured process output for target model output training data, and by:
(a) setting an estimated dominant time constant associated with each input training data;
(b) for each input training data, initializing a filter structure based upon the input training data's estimated dominant time constant;
(c) for each target device output training data, determining a corresponding target analyzer output training data (d) mapping input training data through the delay units and the filter structure, resulting in a plurality of training filter state vectors, creating a linear or non-linear analyzer trained to map training states to the corresponding target analyzer outputs, where said analyzer has a set of internal states which are linearly dependent on a current filter state vector, resulting in a linear multi-input, multi-output state space system mapping input training data to said internal analyzer states, where said linear multi- input, multi-output state space system consists of a block structure of single-input multi-output sub-systems, one for each filter unit; (e) applying a method of order reduction to each single-input, multi-output state space system, resulting in a set of reduced order single-input, multi- output sub-systems; and
(f) configuring each filter unit using the corresponding reduced order single- input, multi-output sub-systems.
34. The method of claim 33 wherein the non-linear analyzer step employs a neural network and the internal states are the summation values of the neural network hidden layer.
35. The method of claim 33 wherein the linear or non-linear analyzer employs a Partial Least Squares model and the internal states are latent variable scores.
36. The method of Claim 32 or Claim 33 wherein the step of balancing employs a Hankel model reduction procedure.
37. The method of Claim 32 or Claim 33 wherein the step of balancing employs an Internally Balanced model reduction procedure.
38. The method of Claim 32 or Claim 33 wherein the step of balancing employs an iterative application of any of a variety of model reduction procedures including a Hankel model reduction procedure and an Internally Balanced model reduction procedure.
39. The method of Claim 27 wherein, in a predictive mode, the step of receiving includes receiving a contiguous stream of control data inputs from .an external system, said data inputs representing measurements from the modeled process; said receiving of data inputs occurs once per base sample; and the steps (a) through (e) are performed once per base sample.
40. The method of Claim 39 further comprising the step of passing the model output to an external system for analysis, monitoring, optimization or control of the modeled process.
41. The method of Claim 27 wherein, in an horizon mode, steps (a) through (e) are iterated multiple times wherein, at each iteration, the filtered states and other state information are stored for later use.
42. The method of Claim 41 further comprising the step of passing the model output at each iteration to an external device or method for analysis, monitoring, optimization or control of the modeled process.
43. The method of Claim 27 wherein, in a reverse horizon mode , steps (a) through (e) are iterated multiple times in reverse order, wherein, at each iteration the steps employ stored information from the method of Claim 41; and the result of this method is to produce sensitivities of a predictive model output for a specified iteration with respect to changes in the predictive mode received data at each previous iteration, where said specified iteration is provided to the method by an external system.
44. The method of Claim 43 further comprising the step of passing the calculated sensitivities at each iteration to an external device or method for analysis, monitoring, optimization or control of the modeled process.
45. The method of Claim 32 or Claim 33 wherein the step of initializing the filter structure uses a Laguerre expansion.
46. The method of Claim 27 wherein the plurality of filter units comprise: first and/or second order subfilters.
47. The method of Claim 27 wherein the non-linear analyzer comprises: a neural network.
48. The method of Claim 27 wherein the non-linear analyzer comprises: a linear or non-linear Partial Least Squares model.
49. The method of Claim 27 wherein the non-linear analyzer comprises: a hybrid parallel combination of a linear model with a non-linear model.
PCT/US1998/020295 1997-10-01 1998-09-29 Non-linear dynamic predictive device WO1999017175A1 (en)

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
US6063897P 1997-10-01 1997-10-01
US60/060,638 1997-10-01
US09/160,128 US6453308B1 (en) 1997-10-01 1998-09-24 Non-linear dynamic predictive device
US09/160,128 1998-09-24

Publications (1)

Publication Number Publication Date
WO1999017175A1 true WO1999017175A1 (en) 1999-04-08

Family

ID=26740150

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/US1998/020295 WO1999017175A1 (en) 1997-10-01 1998-09-29 Non-linear dynamic predictive device

Country Status (2)

Country Link
US (2) US6453308B1 (en)
WO (1) WO1999017175A1 (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2003036396A1 (en) * 2001-10-25 2003-05-01 Aspen Technology, Inc. Non-linear dynamic predictive device
US7065511B2 (en) 1997-10-01 2006-06-20 Aspen Technology, Inc. Non-linear dynamic predictive device

Families Citing this family (64)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1014240A1 (en) * 1998-12-17 2000-06-28 Siemens Aktiengesellschaft A system of case-based reasoning for sensor prediction in a technical process, especially in a cement kiln, method and apparatus therefore
US7085688B1 (en) * 1999-10-22 2006-08-01 Shizuo Sumida Non-linear characteristic reproducing apparatus and non-linear characteristic reproducing program storage medium
US6826521B1 (en) * 2000-04-06 2004-11-30 Abb Automation Inc. System and methodology and adaptive, linear model predictive control based on rigorous, nonlinear process model
WO2002003152A2 (en) * 2000-06-29 2002-01-10 Aspen Technology, Inc. Computer method and apparatus for constraining a non-linear approximator of an empirical process
JP3806637B2 (en) * 2001-10-30 2006-08-09 三菱電機株式会社 Semiconductor pressure detector
AU2002357828A1 (en) * 2001-12-14 2003-06-30 The Texas A And M University System System for actively controlling distributed applications
US6901300B2 (en) * 2002-02-07 2005-05-31 Fisher-Rosemount Systems, Inc.. Adaptation of advanced process control blocks in response to variable process delay
US7089534B2 (en) * 2002-05-01 2006-08-08 International Business Machines Corporation Model based test generation for validation of parallel and concurrent software
SE522691C3 (en) * 2002-06-12 2004-04-07 Abb Ab Dynamic on-line optimization of production processes
EP1404046B1 (en) * 2002-09-26 2006-11-08 Lucent Technologies Inc. Equalizer with two space-time filters and with a selector for choosing the filter with the best symbol estimate
US20050187643A1 (en) * 2004-02-19 2005-08-25 Pavilion Technologies, Inc. Parametric universal nonlinear dynamics approximator and use
US7529703B2 (en) * 2003-11-18 2009-05-05 Citigroup Global Markets, Inc. Method and system for artificial neural networks to predict price movements in the financial markets
US7451003B2 (en) * 2004-03-04 2008-11-11 Falconeer Technologies Llc Method and system of monitoring, sensor validation and predictive fault analysis
US20050228150A1 (en) * 2004-04-13 2005-10-13 Union Carbide Chemicals & Plastics Technology Corporation Use of instantaneous split to improve reactor control
US7437014B2 (en) * 2004-04-21 2008-10-14 The Boeing Company Methods and systems for model reduction and system identification of dynamic systems with multiple inputs
US7460916B2 (en) * 2004-10-19 2008-12-02 Optichron, Inc. Nonlinear system observation and control
US7467614B2 (en) 2004-12-29 2008-12-23 Honeywell International Inc. Pedal position and/or pedal change rate for use in control of an engine
US8380773B2 (en) * 2005-02-18 2013-02-19 Netlogic Microsystems, Inc. System and method for adaptive nonlinear filtering
US7444193B2 (en) 2005-06-15 2008-10-28 Cutler Technology Corporation San Antonio Texas (Us) On-line dynamic advisor from MPC models
US7389773B2 (en) 2005-08-18 2008-06-24 Honeywell International Inc. Emissions sensors for fuel control in engines
US7451004B2 (en) * 2005-09-30 2008-11-11 Fisher-Rosemount Systems, Inc. On-line adaptive model predictive control in a process control system
US7444272B2 (en) * 2005-10-13 2008-10-28 International Business Machines Corporation Self-modulation in a model-based automated management framework
US7949417B2 (en) * 2006-09-22 2011-05-24 Exxonmobil Research And Engineering Company Model predictive controller solution analysis process
US8671065B2 (en) * 2007-05-29 2014-03-11 Commonwealth Scientific And Industrial Research Organisation Methods and apparatus of monitoring an evolving system using selected functional nest
US8060290B2 (en) 2008-07-17 2011-11-15 Honeywell International Inc. Configurable automotive controller
FR2939928B1 (en) * 2008-12-15 2012-08-03 Snecma STANDARDIZATION OF DATA USED FOR MONITORING AN AIRCRAFT ENGINE
CA2754690C (en) * 2009-03-11 2017-01-24 Exxonmobil Upstream Research Company Gradient-based workflows for conditioning of process-based geologic models
EP2406750B1 (en) 2009-03-11 2020-04-01 Exxonmobil Upstream Research Company Adjoint-based conditioning of process-based geologic models
US8620461B2 (en) 2009-09-24 2013-12-31 Honeywell International, Inc. Method and system for updating tuning parameters of a controller
US8538901B2 (en) 2010-02-05 2013-09-17 Toyota Motor Engineering & Manufacturing North America, Inc. Method for approximation of optimal control for nonlinear discrete time systems
US8504175B2 (en) * 2010-06-02 2013-08-06 Honeywell International Inc. Using model predictive control to optimize variable trajectories and system control
US9335042B2 (en) * 2010-08-16 2016-05-10 Emerson Process Management Power & Water Solutions, Inc. Steam temperature control using dynamic matrix control
US9217565B2 (en) * 2010-08-16 2015-12-22 Emerson Process Management Power & Water Solutions, Inc. Dynamic matrix control of steam temperature with prevention of saturated steam entry into superheater
US9447963B2 (en) 2010-08-16 2016-09-20 Emerson Process Management Power & Water Solutions, Inc. Dynamic tuning of dynamic matrix control of steam temperature
US8489632B1 (en) * 2011-06-28 2013-07-16 Google Inc. Predictive model training management
US9170572B2 (en) * 2011-07-06 2015-10-27 Honeywell International Inc. Dynamic model generation for implementing hybrid linear/non-linear controller
US8649884B2 (en) * 2011-07-27 2014-02-11 Honeywell International Inc. Integrated linear/non-linear hybrid process controller
US9677493B2 (en) 2011-09-19 2017-06-13 Honeywell Spol, S.R.O. Coordinated engine and emissions control system
US9163828B2 (en) 2011-10-31 2015-10-20 Emerson Process Management Power & Water Solutions, Inc. Model-based load demand control
US9650934B2 (en) 2011-11-04 2017-05-16 Honeywell spol.s.r.o. Engine and aftertreatment optimization system
US20130111905A1 (en) 2011-11-04 2013-05-09 Honeywell Spol. S.R.O. Integrated optimization and control of an engine and aftertreatment system
US9519858B2 (en) * 2013-02-10 2016-12-13 Microsoft Technology Licensing, Llc Feature-augmented neural networks and applications of same
US20140278285A1 (en) * 2013-03-13 2014-09-18 Vasilis Z. Marmarelis Method for Quantitative Diagnosis of Cerebrovascular, Neurovascular and Neurodegenerative Diseases via Computation of a CO2 Vasomotor Reactivity Index based on a Nonlinear Predictive Model
CN103941589B (en) * 2014-04-24 2016-08-24 中国科学院自动化研究所 A kind of nonlinear model predictive control method of piezo actuator
KR102274069B1 (en) * 2014-10-30 2021-07-06 삼성에스디에스 주식회사 Apparatus and method for generating prediction model
EP3051367B1 (en) 2015-01-28 2020-11-25 Honeywell spol s.r.o. An approach and system for handling constraints for measured disturbances with uncertain preview
EP3056706A1 (en) 2015-02-16 2016-08-17 Honeywell International Inc. An approach for aftertreatment system modeling and model identification
EP3091212A1 (en) 2015-05-06 2016-11-09 Honeywell International Inc. An identification approach for internal combustion engine mean value models
EP3125052B1 (en) 2015-07-31 2020-09-02 Garrett Transportation I Inc. Quadratic program solver for mpc using variable ordering
US10272779B2 (en) 2015-08-05 2019-04-30 Garrett Transportation I Inc. System and approach for dynamic vehicle speed optimization
JP6443311B2 (en) * 2015-11-30 2018-12-26 オムロン株式会社 Control device, control program, and recording medium
US10415492B2 (en) 2016-01-29 2019-09-17 Garrett Transportation I Inc. Engine system with inferential sensor
US10124750B2 (en) 2016-04-26 2018-11-13 Honeywell International Inc. Vehicle security module system
US10036338B2 (en) 2016-04-26 2018-07-31 Honeywell International Inc. Condition-based powertrain control system
WO2018101918A1 (en) 2016-11-29 2018-06-07 Honeywell International Inc. An inferential flow sensor
US9899013B1 (en) * 2017-03-28 2018-02-20 EarthQuaker Devices, LLC Forward and reverse delay effects pedal
US11057213B2 (en) 2017-10-13 2021-07-06 Garrett Transportation I, Inc. Authentication system for electronic control unit on a bus
US11387790B2 (en) 2019-01-24 2022-07-12 Analog Devices International Unlimited Company Power semiconductor device with charge trapping compensation
US10985951B2 (en) 2019-03-15 2021-04-20 The Research Foundation for the State University Integrating Volterra series model and deep neural networks to equalize nonlinear power amplifiers
RU2720219C1 (en) * 2019-06-07 2020-04-28 Федеральное государственное бюджетное образовательное учреждение высшего образования "Казанский государственный энергетический университет" (ФГБОУ ВО "КГЭУ") Adaptive digital predictive device
RU2720218C1 (en) * 2019-06-07 2020-04-28 Федеральное государственное бюджетное образовательное учреждение высшего образования "Казанский государственный энергетический университет" (ФГБОУ ВО "КГЭУ") Digital predictive device
US11533070B2 (en) * 2019-12-23 2022-12-20 Analog Devices International Unlimited Company Systems and methods of compensating for narrowband distortion in power semiconductor devices
US11563409B2 (en) 2020-10-26 2023-01-24 Analog Devices International Unlimited Company Configurable non-linear filter for digital pre-distortion
WO2022169586A1 (en) * 2021-02-02 2022-08-11 Flex Logix Technologies, Inc. Mac processing pipeline having activation circuitry, and methods of operating same

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO1997028669A1 (en) * 1996-01-31 1997-08-07 Asm America, Inc. Model-based predictive control of thermal processing
US5659667A (en) * 1995-01-17 1997-08-19 The Regents Of The University Of California Office Of Technology Transfer Adaptive model predictive process control using neural networks

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6823244B2 (en) * 1995-06-07 2004-11-23 Automotive Technologies International, Inc. Vehicle part control system including electronic sensors
US5535135A (en) * 1993-08-24 1996-07-09 Motorola, Inc. State estimator based exhaust gas chemistry measurement system and method
US5877954A (en) * 1996-05-03 1999-03-02 Aspen Technology, Inc. Hybrid linear-neural network process control
US5992383A (en) * 1996-05-28 1999-11-30 U.S. Philips Corporation Control unit having a disturbance predictor, a system controlled by such a control unit, an electrical actuator controlled by such a control unit, and throttle device provided with such an actuator
US6453308B1 (en) 1997-10-01 2002-09-17 Aspen Technology, Inc. Non-linear dynamic predictive device
US6493689B2 (en) * 2000-12-29 2002-12-10 General Dynamics Advanced Technology Systems, Inc. Neural net controller for noise and vibration reduction

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5659667A (en) * 1995-01-17 1997-08-19 The Regents Of The University Of California Office Of Technology Transfer Adaptive model predictive process control using neural networks
WO1997028669A1 (en) * 1996-01-31 1997-08-07 Asm America, Inc. Model-based predictive control of thermal processing

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
B.GIBBS ET AL: "NONLINEAR MODEL PREDICTIVE CONTROL FOR FOSSIL POWER PLANTS", PROCEEDINGS OF THE 1992 AMERICAN CONTROL CONFERENCE, vol. 4, 24 June 1992 (1992-06-24), USA, pages 3091 - 3098, XP002091769 *
M.KATEBI ET AL: "PREDICTIVE CONTROL DESIGN FOR LARGE SCALE SYSTEMS", INTEGRATED SYSTEMS ENGINEERING.A POSTPRINT VOLUME FROM THE IFAC CONFERENCE, 27 September 1994 (1994-09-27), UK, pages 47 - 52, XP002091768 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7065511B2 (en) 1997-10-01 2006-06-20 Aspen Technology, Inc. Non-linear dynamic predictive device
WO2003036396A1 (en) * 2001-10-25 2003-05-01 Aspen Technology, Inc. Non-linear dynamic predictive device

Also Published As

Publication number Publication date
US6453308B1 (en) 2002-09-17
US7065511B2 (en) 2006-06-20
US20020178133A1 (en) 2002-11-28

Similar Documents

Publication Publication Date Title
US7065511B2 (en) Non-linear dynamic predictive device
Hou et al. From model-based control to data-driven control: Survey, classification and perspective
Maciejowski et al. Predictive control
US8046090B2 (en) Apparatus and method for automated closed-loop identification of an industrial process in a process control system
Andrieu et al. Particle methods for change detection, system identification, and control
US8055358B2 (en) Multi-objective predictive process optimization with concurrent process simulation
US9494925B2 (en) Hybrid control system
US20070078529A1 (en) On-line adaptive model predictive control in a process control system
US20050187643A1 (en) Parametric universal nonlinear dynamics approximator and use
CA2678820A1 (en) Two-stage model predictive control technique
Wang A Tutorial on Model Predictive Control: Using a Linear Velocity‐Form Model
WO2003036396A1 (en) Non-linear dynamic predictive device
Faraji-Niri et al. Stochastic stability and stabilization of semi-Markov jump linear systems with uncertain transition rates
Giulioni Stochastic model predictive control with application to distributed control systems
Emhemed et al. Model predictive control: A summary of industrial challenges and tuning techniques
Miranda et al. Hierarchical approach for ℋ∞ robust control design: S/KS mixed sensitivity with genetic algorithm
Reinelt H∞ loop shaping for systems with hard bounds
Van Willigenburg et al. UDU factored discrete-time Lyapunov recursions solve optimal reduced-order LQG problems
Witczak et al. Observers and genetic programming in the identification and fault diagnosis of non-linear dynamic systems
Zhu et al. Analysis and tuning of a three-term DMC
Huzmezan et al. Multivariable Laguerre-based indirect adaptive predictive control a reliable practical solution for process control
Rahim Detection and elimination of insignificant interacting subsystems in MIMO closed-loop systems using the least mean square-based partial correlation algorithm
Rizvi Reinforcement Learning for Model-Free Output Feedback Optimal Control
Dong et al. Multimode Diagnosis for Switched Affine Systems
Duddeti Approximation of Fractional-Order Systems Using Balanced Truncation with Assured Steady-State Gain

Legal Events

Date Code Title Description
AK Designated states

Kind code of ref document: A1

Designated state(s): CA JP

AL Designated countries for regional patents

Kind code of ref document: A1

Designated state(s): AT BE CH CY DE DK ES FI FR GB GR IE IT LU MC NL PT SE

121 Ep: the epo has been informed by wipo that ep was designated in this application
DFPE Request for preliminary examination filed prior to expiration of 19th month from priority date (pct application filed before 20040101)
122 Ep: pct application non-entry in european phase
NENP Non-entry into the national phase

Ref country code: CA