A NOVEL ALGORITHM FORMODEL PLANT MISMATCHDETECTION FOR MODELPREDICTIVE CONTROLLERSbyYiting TsaiB.A.Sc., The University of British Columbia, 2013A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Chemical and Biological Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)April 2016c© Yiting Tsai 2016AbstractFor Model Predictive Controlled (MPC) plants, the quality of the plant model determinesthe quality of performance of the controller. Model Plant Mismatch (MPM), the discrep-ancies between the plant model and actual plant transfer matrix, can both improve ordegrade controller “performance,” depending on the context in which“performance” is de-fined. Instead of using simple, “yes-no” type of performance metrics to diagnose whetherMPM is present, this thesis achieves the further goal of pinpointing the exact plant inputscausing MPM. The detection method consists of a two-step identification procedure. Thefirst experiment identifies the MPM-affected rows in the plant matrix. The second exper-iment pinpoints the exact inputs causing said MPM, and these MPM-affected elementsare further ascertained using a hypothesis test. Finally, using input design, an estimatedoptimal excitation sequence is generated based on the available closed-loop data, whichhelps the engineer determine whether the original excitation was sufficient for the afore-mentioned sys-ID experiments. The most important underlying assumptions are the lineartime-invariance of the plant and noise dynamics. The proposed algorithm is exercised onartificial 3x3 and 5x5 plants, then on real data from an industrial lime slaker suffering fromsparse MPM, to demonstrate its ability of accurately pinpointing MPM-affected elements.iiPrefaceThis thesis was written to fulfill the requirements of the degree of Master of Applied Science,in the Faculty of Graduate and Postdoctoral Studies, in the Department of Chemical andBiological Engineering at the University of British Columbia. The contents describe a novel,System Identification-based approach to detect both the presence of MPM and the exacttransfer elements affected by MPM, using routine closed-loop data. This saves significanttime and production downtime, because the amount of disruption introduced upon theprocess is minimized.The first paper on MPM detection, under the special case of linear time-invariance andwhite noise disturbance, was submitted to the 2015 International Symposium on AdvancedControl of Chemical Processes (ADCHEM). This paper explored the basic System Identi-fication methods to detect both the existence and locations of MPM, but did not providea mathematical way of justifying the locations, or whether the setpoint was persistentlyexciting. This thesis addresses the issue of MPM location pinpointing using hypothesistesting, and determines an optimal excitation signal using input design.I, Yiting Tsai, was the main contributor for the theoretical research and developmentof this novel MPM detection algorithm. Bhushan Gopaluni, my supervisor, provided in-valuable guidance in most areas of the research. Finally, Spartan Controls provided bothfunding and industrial expertise in order to keep the development relevant to the industri-ally available process control monitoring and analysis tools.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Model-Predictive Control (MPC) . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Identification and Modelling of LTI Systems . . . . . . . . . . . . . . . . . . 71.3 Delay Estimation of LTI Systems . . . . . . . . . . . . . . . . . . . . . . . . 142 Model Plant Mismatch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.1 MPM Detection based on Performance Assessment . . . . . . . . . . . . . . 21ivTABLE OF CONTENTS2.2 MPM Detection based on Correlations . . . . . . . . . . . . . . . . . . . . . 242.3 MPM Detection based on Partial Correlations . . . . . . . . . . . . . . . . . 273 The Proposed MPM Detection Algorithm . . . . . . . . . . . . . . . . . . 303.1 Detecting the Existence of Row MPM: The S∆ Test . . . . . . . . . . . . . 323.2 Pinpointing the Specific Locations of MPM . . . . . . . . . . . . . . . . . . 373.3 Time Delay Estimation for LTI Models . . . . . . . . . . . . . . . . . . . . . 393.4 Hypothesis Testing for “Significant” Model Coefficients . . . . . . . . . . . . 413.5 Input Design for the Excitation Signal . . . . . . . . . . . . . . . . . . . . . 453.6 Computational Savings of the Proposed Method . . . . . . . . . . . . . . . 534 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564.1 3x3 Plant: Case of No MPM . . . . . . . . . . . . . . . . . . . . . . . . . . 564.2 3x3 Plant: Case of Sparse MPM . . . . . . . . . . . . . . . . . . . . . . . . 584.3 5x5 Plant: Case of Sparse MPM . . . . . . . . . . . . . . . . . . . . . . . . 625 Application of MPM Detection Algorithm on Industrial Plant Data . . 696 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 786.1 Summary of Present Contributions . . . . . . . . . . . . . . . . . . . . . . . 786.2 Areas for Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82vList of Tables1.1 List of commonly used linear models . . . . . . . . . . . . . . . . . . . . . . 113.1 Computational Requirement Comparison . . . . . . . . . . . . . . . . . . . 533.2 Computational Requirement Comparison . . . . . . . . . . . . . . . . . . . 544.1 BJ model between d1 and u . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.2 BJ model between d2 and u . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.3 Bi,j,max coefficients for 3x3 sparse MPM case . . . . . . . . . . . . . . . . . 614.4 Hypothesis testing for coefficients for 3x3 sparse MPM case (α = 1%) . . . 614.5 5x5 plant: BJ model between d1 and u . . . . . . . . . . . . . . . . . . . . . 654.6 5x5 plant: BJ model between d2 and u . . . . . . . . . . . . . . . . . . . . . 654.7 5x5 plant: BJ model between d4 and u . . . . . . . . . . . . . . . . . . . . . 664.8 Bi,j,max coefficients for 5x5 sparse MPM case . . . . . . . . . . . . . . . . . 674.9 Hypothesis testing for coefficients for 5x5 sparse MPM case (α = 1%) . . . 675.1 Bi,j,max coefficients for industrial lime slaker . . . . . . . . . . . . . . . . . . 755.2 Standard deviations for Bi,j,max coefficients for industrial lime slaker . . . . 765.3 Hypothesis testing for coefficients for industrial lime slaker . . . . . . . . . . 76viList of Figures1.1 Setpoint improvement due to MPC . . . . . . . . . . . . . . . . . . . . . . . 31.2 The Receding Horizon concept . . . . . . . . . . . . . . . . . . . . . . . . . 41.3 First-Order CSTR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.4 A train of CSTRs with both series and parallel connections . . . . . . . . . 91.5 A black box disregarding all physical relationships between input(s) andoutput(s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.6 Choosing model parameter values based on the minimization of (a functionof) simulation errors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.1 Basic interpretation of MPM . . . . . . . . . . . . . . . . . . . . . . . . . . 202.2 SISO IMC structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.3 SISO IMC structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.1 Internal Model Controller (IMC) Structure . . . . . . . . . . . . . . . . . . 323.2 Combining S∆ and ∆ methods to pinpoint exact locations of MPM . . . . . 393.3 Probability distribution and confidence interval of Bi,j,max . . . . . . . . . . 443.4 Typical excitation signals used in industry . . . . . . . . . . . . . . . . . . . 454.1 Overall S∆(ω) plot for the no-MPM case . . . . . . . . . . . . . . . . . . . 584.2 Overall and row S∆(ω) plots for the sparse MPM case . . . . . . . . . . . . 594.3 Optimal spectrum and setpoint signals for MPM Detection . . . . . . . . . 62viiLIST OF FIGURES4.4 Overall and row S∆(ω) plots for a 5x5 plant suffering from sparse MPM. . . 644.5 Optimal spectrum and setpoint signals for MPM Detection . . . . . . . . . 685.1 Kraft process for the delignification of wood chips . . . . . . . . . . . . . . . 705.2 Lime slaker for the conversion of lime to calcium hydroxide. . . . . . . . . . 715.3 Normalized process variables for the lime slaker . . . . . . . . . . . . . . . . 735.4 Normalized slaker temperature rise y1 (red) and slaker temperature risesetpoint r1 (black) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 745.5 ||S∆||∞ values for the lime slaker . . . . . . . . . . . . . . . . . . . . . . . . 755.6 Optimal spectrum and setpoint signals for MPM Detection . . . . . . . . . 77viiiNomenclatureARX Auto-Regressive, Extra InputsARMA Auto-Regressive, Moving AverageARMAX Auto-Regressive, Extra Inputs, Moving AverageBJ Box-JenkinsC Controller transfer matrixCV Control Variable∆ Additive Mismatch Matrixd Simulation error signalse Zero-mean white noiseDV Disturbance Variable Controller input signalsFIR Finite Impulse ResponseFOPTD First-Order Plus Time Delay Transfer Functionh Finite Impulse Response coefficientsH Finite Impulse Response transfer functionH∞ Hardy-space infinite norm; maximum singular valueLn nth Lebesgue-space normLTI Linear Time-InvariantMA Moving AverageixNomenclatureMIMO Multiple-Input, Multiple-OutputMPM Model Plant MismatchMV Manipulated VariableN(0, α) Gaussian distribution with zero mean and variance αOE Output ErrorP Plant transfer matrixPˆ Plant model transfer matrixPE simulation errorPEM simulation error MinimizationPID Proportional-Integral-Derivative ControlΦr Setpoint power spectrumΦu Input power spectrumr Controller setpoint signalsS∆ Model matrix between setpoints and control errorsSISO Single-Input, Single-OutputSNR Signal-to-Noise RatioSys− ID System Identificationu Plant input signalsy Plant output signalsv Disturbance signalsyˆ Plant model output signalsxAcknowledgementsI would like to extend profound gratitude towards Spartan Controls, for providing prac-tical insight and expertise with cutting-edge process control software and methods usedin industry, NSERC for its generous financial support for this project, and my supervisorBhushan Gopaluni for his invaluable guidance during my research. Finally, I offer specialthanks to my parents, especially my father, for their unyielding support throughout theyears of my studies.xiDedicationTo my parentsxiiChapter 1Introduction1.1 Model-Predictive Control (MPC)Before the rise of MPC, Proportional-Integral-Derivative (PID) control was widely used,since it allowed controllers to respond to spontaneous events that occurred during plant op-eration. PID control was invaluable, because without it, operators would have to respondto abnormal events by changing process variables in an ad hoc manner. For example,consider a plant consisting of a train of reactors, and an abrupt change of process con-ditions occurs in an upstream reactor. Operators must respond quickly by changing theprocess variables in the downstream reactors, in order to maintain an acceptable level ofproduction. However, the control actions suggested by PID feedback control are often sub-optimal, because the control strategy does not account for the existence of various processconstraints, specifically ones that accounted for physical limitations, safety hazards, andeconomic considerations. These constraints further limit the ranges of process variablesallowed; if MPC is used, optimal control actions that remained within the constraint(s)can be easily determined. The most common consequences of sub-optimal control action(due to the lack of MPC) are controller performance degradation, and, in some severecases, plant instabilities [13].In traditional (non-MPC-based) automatic control, setpoints and various constraintsexist for relevant process variables [17]. Setpoints are the desired values of Controlled11.1. Model-Predictive Control (MPC)Variables (CVs), which must be maintained around certain values (the most commonexamples being the level of a fluid in a tank or the temperature of a reactor). These valuesare often determined and justified by the maximization of profit obtained from operating atsuch values, and may change dynamically over time. Setpoints are maintained around theirdesired values by changing the Manipulated Variables (MVs), whose values can be altereddirectly by plant engineers through mechanical means such as valves. Constraints on theprocess typically consist of physical, safety, and economic considerations, as mentionedpreviously. Typical PID control is capable of maintaining CVs around their setpoints.However, one of the following two scenarios often occur as a result. In the first case, theconstraints are not violated, but the controller suggests a mode of operation far below theconstraints that hinders the productivity of the plant. In the second case, the controller issuggesting a mode of operation that produces adequately, but the hard physical constraintsand/or safety constraints are violated. MPC is superior over traditional control strategiesdue to its ability to predict future CV trajectories, and suggest control modes close toor at the constraints of the plant. The benefits of doing so is illustrated in the followingexample.Consider a chemical process, an exothermic reactor, operating at a temperature T .Obviously, T must be less than some maximum allowable temperature, Tmax, in order toprevent exothermic runaway. Due to changes in MVs, plant operators maintain T at valuesmuch smaller than Tmax to minimize the probability of violating the constraint. Figure 1.1illustrates the hypothetical, bell-shaped probabilistic distributions of temperatures.21.1. Model-Predictive Control (MPC)Figure 1.1: Setpoint improvement due to MPCFor traditional PID control, the blue curve represents the temperature distribution,and Told is the average operational temperature that minimizes the probability of violatingTmax. Although the constraint (Tmax) is obviously not violated, a significant gap betweenTold and Tmax exists. The plant operator would wish to increase the average operationaltemperature, in order to increase productivity. By using MPC, the temperature distribu-tion is narrowed and the average operational temperature is moved to closer to Tnew, avalue closer to Tmax. Therefore, MPC achieves the same goal of minimizing the probabilityof violating Tmax, while simultaneously maximizing the production rate.The theory behind MPC is based on the Receding Horizon concept, similar to all othervariants of Predictive Control. Consider the following signal trajectories of a hypotheticalSingle-Input, Single-Output (SISO) plant:31.1. Model-Predictive Control (MPC)Figure 1.2: The Receding Horizon conceptIn Figure 1.1, the relevant signals include the input u(n) and output y(n), with theindex n representing the sampling interval. Denoting the value k as the current samplinginterval, the current input and output values are u(k) and y(k). The plant setpoints,set(n), vary over time and follow the green trajectory. Due to system dynamics, theoutput y(n) cannot reach set(n) immediately, and will instead follow an ideal trajectorycloser to r(n|k), defined as the Ideal Reference Trajectory evaluated at sampling intervaln = k. Note that the Reference Trajectory is re-defined and re-calculated for every samplinginterval n = k, n = k + 1, n = k + 2, etc, given the trajectories of past inputs and outputs.The span of future output samples that the controller is asked to predict is known as thePrediction Horizon, denoted by Hp. Similarly, the span of control action is known as theInput Horizon, denoted by Hu (meaning the input will change from k to k+Hu, then stayconstant afterwards). Suppose that the calculated output, yˆ(n|k), is to meet the ReferenceTrajectory at certain sampling times, termed the Coincidence Points. In this example,41.1. Model-Predictive Control (MPC)suppose that only one such Coincidence Point exists at n = k + Hp. The goal of MPC isto calculate an optimal input trajectory, uˆ(n|k), such that yˆ(n|k) will meet the referencetrajectory at the Coincidence Point, while maintaining all relevant process variables withintheir respective constraints. After uˆ(n|k) is determined, the controller will then apply onlythe first input, uˆ(k|k), to the actual process. Then, the optimal trajectory from the startingpoint n = k + 1 to the same coincidence point n = k + Hp is re-calculated from scratch,with the prediction window moved to the interval k + 1 to n = k + 1 +Hp. MPC repeatsthe exercise of calculating the optimal input trajectory uˆ, and only applying the first step,until time interval n = k + Hu is reached. After this sampling interval, the input is heldat the last computed optimal value, uˆ(k + Hu|k + Hu). The “retreating” nature of theprediction window gives the algorithm its fitting name, “Receding Horizon.”The optimal trajectory, uˆ(k + i|k), is calculated assuming that the system can beapproximated as Linear, Time-Invariant (LTI ) in a small range about a specified operatingpoint. This is usually true for most electrical and mechanical processes, as well as severalchemical processes that have well-known dynamics derived from first principles. A costfunction at any sampling time k, V (k), can be defined as the sum of squared simulationerrors. The magnitude of this error represents the controller offset, the difference betweenthe reference trajectory and the calculated output. In order to achieve an acceptablequality of control, the controller offset must be minimized. This minimization exercise isusually expressed as the following optimization problem:51.1. Model-Predictive Control (MPC)minimizeuV (k) =Hp∑i=Hw||r(k + i|k)− yˆ(k + i|k)||2 +Hu−1∑i=0|| 4 uˆ(k + i|k)||2subject to Luˆ ≤ uˆ(k + i|k) ≤ UuˆL4uˆ ≤ 4uˆ(k + i|k) ≤ U4uˆLyˆ ≤ yˆ(k + i|k) ≤ UyˆAInEq · yˆ ≤ bInEq(1.1)Hw is the time at which inputs are applied to the actual process, and 4uˆ(k + i|k) =uˆ(k + i|k) − uˆ(k + i − 1|k) is the difference between inputs to be applied at successivesampling times. The constraints include Uuˆ and Uyˆ, the respective upper bounds on theinputs and outputs, Luˆ and Lyˆ, the respective lower bounds on the inputs and outputs, andL4uˆ and ≤ U4uˆ, the respective lower and upper bounds on the rates of change of inputs.Finally, AInEq and bInEq are the inequality matrix and column vector, respectively, for theinequality constraints.If the system is LTI, and the cost function V (k) is strictly convex, then the minimiza-tion problem can be readily solved using existing numerical solvers for convex functions.Moreover, V (k) can be re-written as a quadratic form by re-defining the simulation erroras eˆ(k + i|k) = rˆ(k + i|k)− uˆ(k + i|k):V (k) =Hp∑i=HweT (k + i|k) ·Q · e(k + i|k) + uT (k + i|k) ·R · u(k + i|k) (1.2)where Q and R denote the Weighting Factors multiplied to the simulation error andinputs.61.2. Identification and Modelling of LTI SystemsIn a typical process plant, the following hierarchy ties together the various levels ofcontrol commands from the top-most level of human set-point inputs, to the “front line” ofthe process where the actuators exist. The economic considerations of the plant decide theplant-wide and unit-dependent setpoints, which usually change on an hourly to daily basisdue to various factors (for example, supply and customer demands, fluctuating chemicalprices, etc.). The bottom-most level of control consists of PI and PID controllers, whicheliminate offsets and achieve basic performance tuning (such as response time, overshoot,etc.) These controllers then send the control signals to the actuators of the process, acommon example being those of valves. Between the PID controllers and the set-point op-timization layers is a series of complex logic overrides, decoupling networks, and exceptionhandling [13] which handles all spontaneously-arising conditions that cannot be handledby the simple PID layer. The Linear Programming Layer solves the MPC cost functionminimization problems, then sends the solutions to the MPC layer. The MPC layer gener-ates a control signal, which acts as a “setpoint” for the PID layer, which then generates thefinal control signal to the actuators. In this fashion, the predictive control layer providessuperiority in performance and (in most cases) stability compared to ad hoc fixes.1.2 Identification and Modelling of LTI SystemsIn most chemical processes, the determination of control actions are based on the modelthat describes the plant dynamics. This model must be well-characterized if acceptablecontrol actions are to be made. If the underlying physics of process variables are well-known, then the plant input-output relationships can be constructed by the underlyingmass/energy balance equations, as well as any physical relationships between variables.This type of model is known as a First-Principles Model. For example, consider a simpleContinuously-Stirred Tank Reactor (CSTR), where a single reactant (A) is transformed71.2. Identification and Modelling of LTI Systemsinto a single product (B) via a first-order chemical reaction, with reaction rate constant k:Ak→ B (1.3)Figure 1.3: First-Order CSTRSuppose that one wishes to find the fractional conversion of reactant A (denoted by X)as a function of the inlet volumetric flowrate of A (denoted by F ). In this case, the plantinput is F and the plant output is X. A mole balance on the CSTR yields the followinginput-output relationship between F and X [7]:X =τ · k1 + τ · k (1.4)In Eq. 1.4 X is expressed as an explicit function of F due to the straightforward massbalance. However, if the problem is transformed into a train of CSTRs with the samechemical reaction kinetics, the relationship between F and X is no longer straightforward,81.2. Identification and Modelling of LTI Systemsand may be impossible to obtain if the system is too large:Figure 1.4: A train of CSTRs with both series and parallel connectionsIn this case, an explicit, physical function of X with respect to F may be found aftera tenuous series of mass balances, but this approach is impractical due to the followingreasons. First, the first-principles model is subject to errors due to its complexity. Secondly,the model is difficult to maintain if a change occuurs in the system, such as a change in thevalue of a physical constant, or the configuration of reactors. Finally, the model is likelyto be too computationally expensive for situations where the control actions (based on theoutput values) must be determined quickly. In light of these shortcomings related to thefirst-principles approach, the problem of determining the input-output relation can insteadbe treated as the identification of black box (which lies between the input and output):91.2. Identification and Modelling of LTI SystemsFigure 1.5: A black box disregarding all physical relationships between input(s) and out-put(s).For the black box, the input data u(t) and output data y(t) are considered, but thephysical interactions between them are ignored completely. Instead, the outputs are mod-elled as functions of the inputs, by choosing from a set of pre-determined, Linear TimeInvariant (LTI) models. These models relate the inputs u(t), outputs y(t), and noise e(t)using linear polynomial coefficients, denoted by capital letters A through F [12]:A(q) · y(t) = B(q)F (q)· u(t) + C(q)D(q)· e(t) (1.5)“q” denotes the forward-shift operator in time-domain (i.e. q · u(t) = u(t + 1), andq−1 · u(t) = u(t− 1)). These model structures are based off three important assumptions:Linear Time-Invariance (LTI), causality, and that the distribution of the noise e(t) is zero-mean Gaussian. Linear Time-Invariance implies that, for the input-output mapping u(t)F→y(t):1. Superposition applies, i.e. F [m1 · u1(t) + m2 · u2(t)] = m1 · F [u1(t)] + m2 · F [u2(t)](where m1, m2 are constants)101.2. Identification and Modelling of LTI Systems2. The mapping F does not depend on the time instance at which it is applied; let Dθrepresent a delay function of θ units, such that Dθ[u(t)] = u(t− θ). Time-invarianceimplies that Dθ{F [u(t)]} = F{Dθ[u(t)]}, i.e. Dθ ◦ F = F ◦Dθ.On the other hand, zero-mean “Gaussian noise” implies that the noise sequence e(t)possesses the following properties:1. The sequence mean value is zero, i.e. E[e(t)] = 02. The sequence possesses a constant variance, i.e. V AR[e(t)] = λ23. The sequence is normally distributed, i.e. e(t) ∈ N(0, λ)4. The sequence is uncorrelated with itself at any two different points in time, i.e.E[e(s)e(t)] = 0 ∀s, t ∈ R | s 6= t.The presence of the polynomial coefficients determine the type of model used; thefollowing complete list can be found in [12]:Table 1.1: List of commonly used linear modelsPolynomial(s) Present Model Type AbbreviationB Finite Impulse Response FIRAB Auto-Regressive, Extra Inputs ARXABC Auto-Regressive, Moving Average, Extra Inputs ARMAXAC Auto-Regressive, Moving Average ARMABF Output Error OEBFCD Box-Jenkins BJFor example, the ARX model can be expanded into the following relationship:y(t) + a1 · y(t− 1) + a2 · y(t− 2) + · · ·+ ana · y(t− na) (1.6)= b1 · u(t− 1) + b2 · u(t− 2) + · · ·+ bnb · u(t− nb) + e(t) (1.7)111.2. Identification and Modelling of LTI Systemswhere na is the oldest output sampling interval, nb the oldest input sampling interval,a1 · · · ana the output polynomial coefficients, and b1 · · · bnb the input polynomial coefficients.In this model, the present output y(t) depends on the regressors, which are a combinedsequence of past sampled inputs and outputs, y(t) · · · y(t − na) u(t) · · ·u(t − nb). Theregressors are usually combined in a vector denoted by φ(t). Moreover, y(t) also containsa contribution from the noise sequence e(t). The polynomials coefficients A, B, C, D, E,and F are also represented as parameter vectors, denoted by θ. They are chosen as valuesthat “best-fit” the observed data, u(t) and y(t). The definition of “best-fit” is familiar indata fitting, but has a more precise definition in control.Figure 1.6: Choosing model parameter values based on the minimization of (a function of)simulation errors.Consider a LTI system which consists of N total datapoints, as well as available inputand output data,u(t) and y(t), combined into a single dataset ZN . A candidate linear modelM is selected to describe the true system, whose real dynamics are unknown (denoted bythe ? sign in the black box). First, a set of input data u(t) is passed through the modelM , and the model output, yˆ(t), is easily determined. Then, the same set of input datau(t) is passed through the real system, which produces the set of real outputs y(t). Thesimulation error d(t) is defined as the difference between the model and actual outputs, i.e.d(t) = y(t) − yˆ(t). The goal in system identification is to select model parameters θ such121.2. Identification and Modelling of LTI Systemsthat a specific function of the simulation error, VN (θ, ZN ) (expressed also as a function ofboth the parameters θ and dataset ZN ), is minimized. Obviously, the number of availablefunctions VN is infinite. However, the most commonly selected function is the well-knownsum of squared errors, expressed as:VN (θ, ZN ) =1NN∑t=112d(t)2 (1.8)Using a re-arrangement of equation 1.6, the model output yˆ(t) can be expressed as aproduct of two vectors:yˆ(t) = [a1 · · · ana b1 · · · bnb ] · [−y(t− 1) · · · − y(t− na) u(t− 1) · · ·u(t− nb)]T = θT · φ(t)(1.9), where θ = [a1 · · · ana b1 · · · bnb ]T represents the model parameters to be estimated, andφ(t) = [−y(t− 1) · · · − y(t−na) − u(t− 1) · · · − u(t−nb)]T represents the regressor vectorcontaining the past input and output measurements. The actual output y(t) will contain,in addition to the model output yˆ(t), errors e(t) due to a combination unmeasured noiseand uncharacterized dynamics, i.e. y(t) = θ · φ(t) + e(t). The minimization of the sum ofsquared errors, VN (θ, ZN ), is known as the least-squares problem:minθ1NN∑t=112[y(t)− yˆ(t)]2 = minθ1NN∑t=112[y(t)− θT · φ(t)]2 (1.10)The parameters θ that minimize VN (θ, ZN ) in the least-squares case is the well-known131.3. Delay Estimation of LTI Systemssolution:θˆLSN = argminθ1NN∑t=112[y(t)− θT · φ(t)]2 = [φ(t)]† · y(t) = [φ(t)Tφ(t)]−1 · φ(t)T y(t)(1.11)The “dagger” † operator represents the famous Moore-Penrose Psuedoinverse fora general matrix of any dimension (i.e. for a matrix A of any size, A† = (AT · A)−1 · AT ).If the noise sequence e(t) is zero-mean Gaussian, then the estimated parameter values θˆLSNare said to be unbiased [12]. The definition of unbiased in the context of control meansthat the expected difference between θˆLSN and the true system parameters θ0 approacheszero as the data size becomes infinitely large:limN→∞E[θˆLSN − θ0] = 0 (1.12)In a similar manner, if the estimated parameters θˆLSN converge to true parameters θ0with infinite data size, then the estimate θˆLSN is said to be consistent.1.3 Delay Estimation of LTI SystemsMany engineering systems possess inherent delays between inputs and outputs, mostly dueto physical factors such as transportation lags or mixing. In rarer cases, physical factorssuch as the malfunction of sensors or measurement instruments can also cause processdelays. Regardless of the cause, an accurate estimate of the delay between a given set ofdata inputs u(t) and outputs y(t) is required for a consistent model.141.3. Delay Estimation of LTI SystemsConsider a nu × ny (i.e. nu inputs, ny outputs) system which assumes a FIR (FiniteImpulse Response) model. In other words, actual outputs y(t) can be expressed as a sumof model outputs yˆ(t) = B(q) · u(t) and zero-mean, Gaussian noise e(t):y(t) = B(q) · u(t) + e(t) (1.13)The individual terms can be expanded in matrix form as the following:y1...yny =B1,1 · · · B1,nu.... . ....Bny ,1 · · · Bny ,nu ·u1...unu+e1...eny (1.14)Each FIR polynomial B contains information on how past inputs affect the currentoutput. Assuming that the system is causal, the current output is only affected by inputsfrom the present sampling interval to the N th earlier sampling interval. In other words, y(t)depends on the inputs u(t), u(t−1), · · · , u(t−N). Therefore, following the model structurein Equation 1.13, the B polynomial corresponding to the ith output and jth input can beexpanded as:Bi,j = Bi,j,0 +Bi,j,1 · q−1 + · · ·Bi,j,N · q−N (1.15)The ith output yi can be expressed as:151.3. Delay Estimation of LTI Systemsyi(t) = Bi,1,0 · u1(t) +Bi,1,1 · u1(t− 1) + · · ·+Bi,1,N · u1(t−N)+Bi,2,0 · u2(t) +Bi,2,1 · u2(t− 1) + · · ·+Bi,2,N · u2(t−N)+...+Bi,nu,0 · unu(t) +Bi,nu,1 · unu(t− 1) + · · ·+Bi,nu,N · unu(t−N)Suppose experimental data Y (t) = [y(t) · · · y(t − N)]T is available. The first outputy(t) contains contributions from the inputs u(t), · · · , u(t−N). A similar dependence holdsfor every other output term, down to the last output y(t−N) (which is dependent on theinputs u(t − N), · · · , u(t − 2N)). Taking advantage of this property, a slightly modifiedregressor matrix Φ can be defined as:Φ(t) =u1(t) · · · u1(t−N) · · · unu(t) · · · unu(t−N)....... . .......u1(t−N) · · · u1(t− 2N) · · · unu(t−N) · · · unu(t− 2N) (1.16)In addition to this, each output term is corrupted by a corresponding noise term (i.e.y(t) by e(t), y(t−1) by e(t−1), etc.). Knowing this, a noise vector E(t) = [e(t) · · · e(t−N)]Tcan also be constructed. On the other hand, the parameters can be collected into a singlevector Θ:Θ = [B1,1,0 · · ·B1,1,N B1,2,0 · · ·B1,2,N · · · B1,nu,0 · · ·B1,nu,N ]T (1.17)161.3. Delay Estimation of LTI SystemsFinally, using the previous vector and matrix definitions, the input-output relationshipcan be expressed as:Y (t) = Φ(t) ·Θ + E(t) (1.18)Θ can be estimated by minimizing a function of the absolute simulation error or residualterm |Φ · Θ − Y |, using an optimization procedure. Then, finding the number of leadingzeroes in the parameters Θ is equivalent to finding the time delay between inputs andoutputs, since a zero coefficient implies a non-significant contribution of the correspondinginput to the output. A viable approach for deterministic systems is to use the unconstrained1-norm minimization [4]:minimizeθ||Θ||1The above algorithm has the advantage of quickly converging to the true solution indeterministic, noise-free systems. In other words, the optimal parameter vector Θ∗ willequal the true parameter vector Θ0. However, the optimization “forces” the coefficientsin the solution to either zero or non-zero, due to the heavy weight placements on small-valued residuals [4]. Therefore, the algorithm estimates the delay with failing accuracy asthe Signal-to-Noise Ratio (SNR) increases.One possible counter is to introduce the quadratic constraint ||Φ · θ − Y ||2 ≤ tol. Thischanges the problem to:171.3. Delay Estimation of LTI Systemsminimizeθ||Θ||1subject to||Φ · θ − Y ||2 ≤ tolThe tolerance tol is a user-specified, educated guess on the upper bound of the magni-tude of noise ||E||2. If a sufficiently sparse, true parameter vector Θ0 exists for the systemsuch that Y = Φ ·Θ0 + E, then the optimal solution Θ∗ will be “sufficiently close” to thetrue parameters Θ0 such that ||Θ∗ − Θ0|| ≤ C · tol, where C is a small constant. Themagnitude of tol will, of course, be chosen based on the application and the convergenceof the optimization problem.18Chapter 2Model Plant MismatchIn any industrial plant containing thousands of process variables, the true plant can neverbe modelled to complete accuracy. The discrepancy between the model and the true systemalways exists, due to the following non-exhaustive list of factors:1. Uncharacterized nonlinearities, or shifts in operation from LTI to non-linear regions2. Disturbances unable to be characterized as a function of Gaussian noise (i.e. CD · e(t)where C,D are linear polynomials in q and e(t) ∈ N(0, λ))3. Model fails to keep up with rapidly-changing process dynamics (can be both linearor non-linear)A famous quote from the prominent 19th-century statistician George Box (who estab-lished the Box-Jenkins model structure) captures this phenomenon perfectly:“All models are wrong; some are useful.” —George Box, 1976.The first part of his statement perfectly captures the concept of Model-Plant Mismatch(MPM): since the true plant is never known exactly, a discrepancy always exists betweenthe model and true system. In the context of process control, MPM is defined as thediscrepant transfer function(s) between the model and true system:19Chapter 2. Model Plant MismatchFigure 2.1: Basic interpretation of MPMThe detection of MPM is important because it may lead to the following inexhaustiblelist of controller and/or plant symptoms:1. Plant Instability: The discrepancy between the real plant and model grows largeenough such that the controller action leads to an unstable response, or the realplant itself becomes unstable and the model fails to diagnose the instabilities in time.2. Performance Degradation: In most cases MPM degrades controller “quality” and“performance.” However, the precise definition of these metrics can vary from userto user, and in some cases the presence of MPM can improve “performance.”3. Disruption of plant operation: Previous approaches used to detect MPM and rectifymodelling errors required the shut-down of normal operation for open-loop Sys-IDexperiments. These disruptions cause significant amounts of lost production and aredetrimental from an economic point of view.Traditionally, MPM is not directly detected. Instead, engineers rely on indirect meth-ods such as performance monitoring or correlations between process variables to measurethe existence of significant MPM. Then, open-loop Sys-ID experiments are performed toupdate the obsolete plant model. Despite the obvious disadvantage of economic loss fromthe halting of plant operation, indirect MPM-detection methods can also be extremely202.1. MPM Detection based on Performance Assessmentunreliable in special cases. These reasons justify a shift towards closed-loop Sys-ID ap-proaches, which are minimally-intrusive in terms of plant operation, and can also providesignificant computational savings depending on the structure of the plant matrices.2.1 MPM Detection based on Performance AssessmentA common indirect approach for MPM-detection is the evaluation of controller performance[18]. In chemical industries, this is usually done in the ad hoc manner described in theIntroduction, where human operators monitor a few important process variables and reactaccordingly when a variable has deviated outside of an acceptable threshold. However, fora relatively complex plant (such as an oil refinery), thousands of variables are involved andthis manual approach is simply not feasible. The performance monitoring task, therefore,would have to be achieved by an automated system which employs certain metrics to judgewhether specific process variables have varied unacceptably.These so-called variations are known as failures in unit operations, instrumentation,and general decline in (chemical) production. They can be caused by either drifts in MVsthemselves, or measured and unmeasured disturbances that were neglected in the processmodel. Despite the cause of the performance degradations, [14] has suggested that theycan be categorized into the following:1. Abrupt : As the name implies, a process variable s abruptly “jumps” to an unaccept-able value, and hovers around said value for a significant period of time. A goodexample of this is if s follows a unit step when plotted over time. This is the mosteasily detectable form of degradation, as its result is often realized as a sudden failurein an unit (such as a valve or pump).2. Incipient : The process variable s does not change abruptly, but instead gradually212.1. MPM Detection based on Performance Assessmentmigrates away from its nominal value over time. This degradation is also easy todetect over time, since the result of the variable drift inevitably leads to some sort offailure in an unit operation.3. Intermittent : The process variable s oscillates (significantly) between a range of finitevalues. This is the most difficult form of degradation to detect, because it is oftenuncertain whether the cause is noise, poor calibration of instruments, poor modellingof system dynamics, or actual faults in variables.Despite these classifications, one cannot easily correlate the observed degradations withthe existence of MPM. Although the degradation in performance is a common symptom ofMPM, a MPM → degradation causal relationship cannot be established for all cases. Onesimple counter-example is illsutrated in the following section. If the performance metricis chosen naively, then the presence of MPM could actually lead to an improvement ofperformance.Example 2.1.1. Consider a SISO, plant with the following Internal Model Controller(IMC) structure [17]:Figure 2.2: SISO IMC structure222.1. MPM Detection based on Performance AssessmentThe relevant signals are the setpoint r(t), real output y(t), model output yˆ(t), andtracking error e(t) = r(t) − y(t). Moreover, suppose the Laplace-domain plant transferfunction is P (s) = 2s+1 , and the plant model is Pˆ (s) =1s+1 . Therefore, the additive MPMis ∆ = P (s)− Pˆ (s) = 1s+1 . Also assume that the controller C(s) is designed as the productof a first-order filter f(s) = 2s+1 and the inverse plant transfer function P−1(s) = s+12 , i.e.C(s) = P−1(s) · f(s) = 1. v(t) is simply some arbitrary noise sequence.Finally, suppose the user defines the controller “performance” as the absolute magni-tude of the tracking error, |e(t)| = |r(t)−y(t)|. Evaluating the actual (in presence of MPM)and designed (in absence of MPM) tracking errors, ea(s) and ed(s), respectively, yields:ea(s) =ss+ 2· [r(s)− v(s)] (2.1)ed(s) =ss+ 1· [r(s)− v(s)] (2.2)If the setpoint experiences a step change, and no disturbance exists, i.e. r(s) = 1s andv(s) = 0, then the absolute tracking errors in time-domain evaluate to:|ea(t)| = e−2t (2.3)|ed(t)| = e−t (2.4)Evidently, MPM has improved the “performance” of the controller in this case, becausethe tracking error in the presence of MPM is smaller than that in the absence of MPM.Clearly, MPM and performance degradation are not tied together in a direct causal232.2. MPM Detection based on Correlationsrelationship. In order to use controller performance as a metric for MPM detection, “per-formance” to be carefully defined in order to avoid a pitfall such as the one previouslyillustrated. In most applications, performance-index detection algorithms are designed ina case-by-case fashion for specific types of plants and controllers, but a universal designthat accounts for all possible combinations of plant dynamics and controller structures isimpossible.2.2 MPM Detection based on CorrelationsThe other common indirect approach of MPM-detection is the measurement of correlationsbetween specific process variables. These methods were mostly introduced and exploredby [20], [1], and [6]. Mathematically, the correlation between two random variables v(t)and w(t) is defined as:Corr(v, w) , Cov(v, w)V AR(v)V AR(w), (2.5)where V AR represents the Variance of a variable:V AR(v) , E[v − E(v)] (2.6)and Cov represents the Covariance between two variables:242.2. MPM Detection based on CorrelationsCov(v, w) , E[v − E(v)][w − E(w)] (2.7)E represents the expected value operator, and N the total number of samples. Thisanalysis can be extended to a MIMO system with m inputs and n outputs, where v and wbecome random vectors such that v ∈ Rm and w ∈ Rn. In this case, the correlation matrixis written as:Corr(v, w) =Corrv1,w1 · · · Corrv1,wn.... . ....Corrvm,w1 · · · Corrvm, wn (2.8)where the individual elements Corr(vi, wj) are defined as in Equation 2.5. Obviously,the variables of interest v and w must be selected such that a significant correlation betweenthem indicates a presence of MPM. Specifically, if the element Corr(vi, wj) is non-zero, thenthe corresponding plant transfer element Pˆij) is experiencing MPM. However, this methodalone is susceptible to pitalls as [6] has demonstrated during open-loop MPM detection,where false positives are diagnosed by the presence of correlations when no actual MPMexists. This is illustrated in the following example.Example 2.2.1. Consider, again, a SISO plant following an IMC structure:252.2. MPM Detection based on CorrelationsFigure 2.3: SISO IMC structureIn this case, the relevant covariances and correlations are evaluated between input u(t)and simulation error d(t):Cov(u, d) = E[u− E(u)][d− E(d)] (2.9)Corr(u, d) =Cov(u, d)V AR(u)V AR(d)(2.10)[20] proved mathematically that if no MPM exists, then Corr(u, d) = 0 (i.e. it isimpossible to find a non-zero correlation between u and d). However, the following examplewill demonstrate that the previous relationship is not commutative. Specifically, if one findsthat Corr(u, d) 6= 0, this does not indicate the presence MPM in all cases.Suppose that a SISO plant exists such that its Laplace-domain model is Pˆ = 1s+1 . Alsoassume that no MPM exists, i.e. Pˆ = P (s) = 1s+1 and ∆(s) = P (s) − Pˆ (s) = 0. Thecontroller C(s) is the product of the plant inverse P−1(s) = s + 1 and first-order filter262.3. MPM Detection based on Partial Correlationsf(s) = 1s+1 , i.e. C(s) = P−1(s) · f(s) = 1. Therefore, the simulation error is equal to thedisturbance, i.e. d(t) = y(t)− yˆ(t) = v(t).Using this relationship, the input u(t) can be re-written as:u(t) = C(q) · [r(t)− d(t)] = r(t)− d(t) = r(t)− v(t) (2.11)The previous equation clearly states that u(t) contains a contribution from d(t) = v(t),therefore Corr(u, d) 6= 0 even though no actual MPM exists. Therefore, a false positivedetection of MPM has been diagnosed.The correlations between the two selected variables of interest, Corr(xi, ej), must beevaluated at all lags. This presents a significant computational burden when selectingcorrelations as the method of choice for MPM detection. Moreover, as shown in the previoussimple example, false detections of MPM (indicated by significant positive correlations) canoccur even when no MPM exists.2.3 MPM Detection based on Partial CorrelationsThis extension to standard correlations was explored by [6], in order to eliminate anypossibility of false positive MPM detections. Suppose that two process variables of interestare v(t) and w(t), collected as time series data. If a third variable z(t) affects both v(t) andw(t), then the correlation between v and w can only be defined for a constant z. In otherwords, if z is allowed to vary arbitrarily, then the correlation between v and w would bemeaningless, because it would be unclear whether a true correlation existed between v andw, or if either one or both of v and w were correlated with z. Therefore, the term Partial272.3. MPM Detection based on Partial CorrelationsCorrelation between v and w is expressed as ParCorr(v, w|z), with the understandingthat z is held constant. [6] uses the following simple method to calculate partial correlation:Bv = z† · v = (zT · z)−1 · zT · v (2.12)rv = v − z ·Bv (2.13)Bw = z† · w = (zT · z)−1 · zT · w (2.14)rw = v − z ·Bw (2.15)ParCorr(v, w|z) = Corr(rv, rw) (2.16)The terms rv and rw are known as Correlation Residuals, and the correlationCorr(rv, rw) is evaluated using Equation 2.5. Expanding this concept to the MIMO casewith nu inputs and ny outputs, partial correlations can be computed for vi and wj for alllags d = 0, 1, 2, · · · , df . Defining a new data vector:v˜i(t) =[u1(t) · · ·unu(t), vi−1(t) · · · vnu(t)]T(2.17), the Partial Correlation can be calculated as:ParCorr(vi(t− d), wj(t)|v˜i(t− d)) (2.18)where holding the previous data entries v˜i(t − d) constant removes the dynamic con-tributions from past data. Just as the case with simple correlations, a significant valuein ParCorr(vi(t − d), wj(t)|v˜i(t − d)) indicates that MPM is present in the plant model282.3. MPM Detection based on Partial Correlationselement Pˆi,j . According to [6], this method yields almost no false positives compared tosimple correlations using closed-loop signal for v and w. The method is foolproof exceptfor the case where the disturbance signal resembles white noise after it passes through thefeedback loop (this pitfall is explored in more detail in [6]). Since the dynamics of eachinput at lag d must be removed to avoid future inputs from generating false positive partialcorrelations, the computational effort involved grows astronomically with the dimensionsof the plant. If the number of inputs and outputs are nu and ny, respectively, and thetotal number of lags df , then [6] estimates the required number of operations to be withinthe order of O(nuny · (nudf )3). Evidently, this approach is only practical with small plantsand/or a small number of lags to be considered.29Chapter 3The Proposed MPM DetectionAlgorithmIn light of the various pitfalls and computation burdens of the traditional performance-based and correlation-based methods of MPM detection, a direct, Sys-ID based approachis considered in this thesis. This method is to provide computational savings, for specifictypes of plants encountered in industry, compared to the traditional model re-identificationperformed when MPM is suspected to exist. The specific objectives are to be accomplishedunder closed-loop operation:1. Detect the existence of MPM with a Sys-ID-based metric2. Pinpoint the exact plant inputs causing MPM, and hence the exact model matrixelements suffering from MPMThe important assumptions listed below suggest the scope (in terms of the type ofplants and controllers) under which this proposed method is applicable:1. MPM is believed to be sparse within the model matrix (i.e. few plant elements aresuspected to be affected by MPM)2. Sufficient setpoint excitation is available, usually in the form of step changes orsinusoids (exact definition of “sufficient excitation” is mathematically established in30Chapter 3. The Proposed MPM Detection AlgorithmSection 3.5)3. Plant dynamics and controller behaviour is stable and approximately LTI within therelevant ranges of operation4. Plant disturbances can be approximated as zero-mean, Gaussian-distributed noise5. Sufficient closed-loop plant data is available for analysisTo reinforce the previous point, this algorithm will not be a suitable candidate for plantsin which most of their model elements are suffering from MPM. In those cases, traditionalsystem re-identification is a superior solution. However, plants which are sparsely affectedby MPM should not have to undergo the trouble of full system re-identification. Themethod of MPM detection should, instead, be one that quickly ascertains the existenceof MPM and identifies the exact model elements which require re-identification. In realapplications, the LTI and sufficient data assumptions may “contradict” each other, in thesense that plant dynamics can only be approximated as LTI in short periods of time,typically in the order of hours to a few days for most chemical plants [13]. If data istaken over long periods of time (more than a few days), then the plant dynamics may havetransitioned to a different linear region, or to a completely non-linear region altogether[14]. In light of these considerations, the control engineer is advised to perform the MPMdetection algorithm using LTI segments of plant data. The frequency of MPM detection candepend on qualitative a priori knowledge of the plant, and/or any empirical informationcollected during operation. For example, consider a hypothetical plant whose dynamicschanges significantly every 6 hours. This can serve as a baseline frequency (i.e. once every6 hours) for detection activation. Moreover, the engineer can also monitor the controllerperformance, and activate MPM detection as soon as the controller performance degradessignificantly. This is the easiest approach to employ in most practical applications, but313.1. Detecting the Existence of Row MPM: The S∆ Testthe performance metric must be carefully designed to avoid the sort of pitfall described inSection 2.1. Finally, the minimum length of data used for detection should be one thatcontains at least one significant setpoint excitation (more rigorously defined in Section 3.5).3.1 Detecting the Existence of Row MPM: The S∆ TestThe familiar Internal Model Control (IMC) structure is again considered, since the plantand plant model transfer matrices are of utmost interest in MPC-related applications:Figure 3.1: Internal Model Controller (IMC) StructureThe relevant signals include r(t) (setpoint), (t) (controller input), u(t) (plant Input),y(t) (real output), yˆ(t) (model output), w(t) (white noise), v(t) (white noise through first-order filter), and d(t) (simulation error).323.1. Detecting the Existence of Row MPM: The S∆ TestMoreover, the relevant transfer functions include C(q) (Controller), P (q) (Real Plant),Pˆ (q) (Plant Model), ∆ = P (q)− ˆP (q) (Additive Model Plant Mismatch), and H(q) (First-Order Noise Filter).The noise sequence v(t) can fall under one of 3 categories:1. Zero-mean, Gaussian-distributed (white)2. First-order coloured (white noise passed through a linear, first-order causal filter)3. Unmeasured noise that cannot be modelled as a combination of white noise andlinear, causal filtersIn most industrial cases, v can usually be approximated as Types 1 or 2. In somespecial cases, drifting disturbances (such as the outdoor temperature for various mills) areconsidered unmeasured, un-modelled noise, unless the drifting behaviour can be empiricallymodelled. Nevertheless, under any type of disturbance v, [2] has derived the relationshipsbetween the relevant IMC signals for the general, discrete-time MIMO case, summarizedas follows:Consider, first, the relationship between the controller input (z), setpoint (z) anddisturbance v(z):(z) = [I + ∆(z) · C(z)]−1 · r(z)− [I + ∆(z) · C(z)]−1 · v(z) (3.1)where I represents the appropriately-sized identity matrix. For simplicity, all followingdiscrete-time sampling arguments (z) will now be omitted from the transfer functions andsignals of interest. The common term S∆ = [I+∆ ·C]−1 is called the “Relative Sensitivity”by [2], because it compares the relative magnitudes of controller errors under the presenceand absence of MPM, respectively. If the “controller error” is defined as the difference333.1. Detecting the Existence of Row MPM: The S∆ Testbetween the setpoint r(z) and real output y(z), i.e. e(z) = r(z) − y(z), then the actual(MPM-present) and design (MPM-absent) control errors (ea and ed, respectively) can beconsidered as:ea = r − y = [I − PˆC][I + ∆C]−1(r − v) (3.2)ed = [I − PˆC](r − v) (3.3)The common term [I − PˆC] can be eliminated from the right sides of the previousequations, by multiplying the left sides by its inverse [I− PˆC]−1. Now the filtered versionsof the actual and design controller errors can be expressed as:eaf = [I − PˆC]−1 · ea = [I + ∆C]−1(r − v) (3.4)edf = [I − PˆC]−1 · ed = r − v (3.5)Finally, the two control errors can be divided against one another to form the ratio ofrelative magnitudes, using the L2-norm. If this norm is evaluated in the frequency domainwithin the control-relevant frequency range Ω, it can be expressed mathematically as:maxω∈Ω,||edf (ω)||2 6=0||eaf (ω)||2||edf (ω)||2= ||S4(ω)||∞ (3.6)The physical interpretation of this ratio can be realized from the fact that control erroris a measure of deviation from perfect setpoint tracking. In other words, if e = r − y 6= 0,then the quality of control is undesirable, depending on the magnitude of the difference.On the other hand, if e = r− y ≈ 0, then the quality of control is acceptable. When MPM343.1. Detecting the Existence of Row MPM: The S∆ Testis absent, the actual and designed control errors are equal in magnitude; this is easilyrealized by substituting ∆ = 0 into Eq. 3.4 and Eq. 3.5. On the other hand, when MPMis present, then the actual and designed control error magnitudes are different. In thiscase, the design (MPM-absent) control error becomes a best-case value to compare to, andtherefore S∆ becomes a sensitivity measure (similar to other engineering contexts) wherethe actual behaviour is compared against a best-case scenario. For this reason, [2] termsthis ratio, which is equivalent to the H∞ norm of the S∆ transfer matrix, the “RelativeSensitivity Index”:RSI , ||S∆(ω)||∞, ω ∈ Ω (3.7)Using this definition, the MPM-absent case (i.e. ∆ = 0) implies that S∆ = I, due to toEq. 3.1. Also, ||S∆(ω)||∞ = ||I(ω)||∞ = 1 due to the equality between eaf and edf . On theother hand, the MPM-present case implies that S∆ 6= I, and ||S∆(ω)||∞ 6= 1. The non-unityvalue of ||S∆(ω)||∞ can therefore serve as an indication of MPM presence. More specifically,Eq. 3.4 and Eq. 3.5 imply that if ||S∆(ω)||∞ > 1, then the actual control error magnitude isgreater than the designed control error magnitude, indicating a degraded overall quality ofcontrol. On the other hand, if ||S∆(ω)||∞ < 1, then the actual control error magnitude issmaller than the designed control error magnitude, indicating an improved overall qualityof control. Despite the possible, interesting implications of these two opposite cases oncontrol quality, they are not explored in this thesis, since the main objective is merelyMPM detection.Eq. 3.1 suggests that S∆ can be obtained by performing Sys-ID between signals (t)and r(t) collected from operating data, and using a suitable LTI model (for example, BJ):353.1. Detecting the Existence of Row MPM: The S∆ TestA(q) · (t) = B(q)F (q)· r(t) + C(q)D(q)· e(t) (3.8)where the error sequence e(t) is assumed to be white noise. If low-order models aredesired, r must be uncorrelated to the disturbance v. On the other hand, if r and v arecorrelated, then a high-order model can be used to obtain an unbiased estimate of S∆ [12].From the time-domain model, a frequency-domain plot of S∆(ω), ω ∈ Ω can be obtained.The RSI represents the H∞ norm, or peak gain on a SISO Bode plot. For MIMO systems,the RSI represents the maximum singular value on a singular value plot.Since each controller input i, i ∈ 1, 2, · · · , ny is independent of all others j , j ∈1, 2, · · · , ny, j 6= i, the user can first identify the entire S∆ matrix, then use the indi-vidual row elements to compute ||S∆(ω)||∞ for each row. This will indicate whether MPMexists in each row of the transfer matrix. In many MPM-absent situations, ||S∆(ω)||∞ willnot exactly equal 1, but should take on values reasonably close to 1. This is due to thequality of the model for the noise term e(t), which may be poor if the disturbance is notLTI.An important note to mention is that the setpoint signal(s) used for the S∆ identi-fication experiment should be sufficiently exciting. Theoretically, in order to identify anunbiased S∆ model of any order ns, the setpoint(s) should be a zero-mean, Gaussian-distributed sequence (with a specified variance) containing all frequencies between 0 tofN (fN being the Nyquist Frequency of the signal, equal to half its sampling rate) [12].Real plants, however, cannot produce excitation signals containing an infinite number offrequencies. Therefore, a good rule of thumb is to use an excitation signal containing ≥ nsnon-zero frequencies, in order to identify an unbiased S∆ model of order ns. If the oppositeis done in practice, i.e. an excitation signal containing < ns non-zero frequencies is used363.2. Pinpointing the Specific Locations of MPMto identify a S∆ model of order ns, the model may become biased. Therefore a trade-offexists between the excitation frequency content and the robustness and reliability of theestimation algorithm, and the control engineer must determine the proper compromisebased on the application at hand.Although S∆ is a powerful metric which can detect the existence of MPM-affectedrows in a MIMO plant matrix, it alone will not suffice if the user wishes to pinpoint theexact MPM-affected elements. This second objective will be accomplished by the methoddescribed in the following section.3.2 Pinpointing the Specific Locations of MPMHaving detected the existence of MPM-affected rows in a general MIMO plant, the nextlogical step is to correctly pinpoint the elements within the MPM-affected rows. Intuitively,the control engineer would wish to conduct Sys-ID experiments to model the MPM matrix,∆. Assume that the MIMO plant has nu inputs and ny outputs, and that the set R containsall the MPM-affected rows. Then, the task of pinpointing is equivalent to identifying thenon-zero elements in ∆Rowi , i ∈ R, since a zero element implies the absence of MPM.To accomplish this, the mathematical relationship between the simulation error d, con-troller error and disturbance v in Figure 3.1, which holds true regardless of the presenceof MPM, can be exploited:d = ∆ · u+ v (3.9)Logically, the transfer function ∆ is no more than a Sys-ID experiment between time-series d and u. [11] pointed out the impossibility of obtaining an unbiased FIR model for∆ using routine, closed-loop data, in the absence of setpoint and/or plant input excitation.373.2. Pinpointing the Specific Locations of MPMHowever, if setpoint excitation is available in the forms of step changes or sinusoids, thenthe Sys-ID experiment for ∆ becomes possible. From the MPM-affected rows singled outby the control engineer using the S∆ test row test, it follows that an existence of MPM inthe jth CV (i.e. j ∈ R) will result in at least one non-zero element in ∆rowj . This is becausedj will contain a contribution from at least one input ui, i ∈ 1, 2, · · · , nu. Conversely, ifno MPM exists in the jth CV, then j /∈ R, and ∆rowj will be entirely zero. One note ofcaution lies on the correlation between d and u. Since the correlation always holds truedue to feedback, a high-order polynomial model must be used for the transfer functionexpression for ∆ [12].For each MPM-affected row, the engineer will use the signals dj , ∀j ∈ R and u1, u2, · · · , unuto estimate unbiased models for each corresponding row of the transfer function ∆(q). Ifthe disturbance v(t) is known to be white noise or first-order filtered white noise, thenhigh-order BJ can be used to obtain an unbiased estimate of ∆(q):A(q) · dj(t) = Bj,1(q)Fj,1(q)· u1(t) + · · ·+ Bj,nu(q)Fj,nu(q)· unu(t) +Cj(q)Dj(q)· vj(t), j ∈ R (3.10)nA, nB, nC , nD, and nF shall be the orders for each polynomial, and each polynomialtakes on the form:M(q) = M1 · q−1 + · · ·+MnM · q−nM (3.11)Each Bj,k polynomial containing significantly non-zero coefficients indicate that dj con-tains a contribution for uk, thus indicating that element ∆j,k contains MPM. Mathemat-ically rigorous definitions of “non-zero” and “sufficient excitation” are developed in detail383.3. Time Delay Estimation for LTI Modelsin Sections 3.4 and 3.5, respectively.The following example illustrates the use of a combination of S∆ and ∆ Sys-ID ex-periments to pinpoint MPM. Suppose the user discovers that ||S∆||∞ = 1 for Row 1, and||S∆||∞ 6= 1 for Rows 2 and 3. Then, Eq. 3.1 implies that Row 1 contains no MPM, andthat Rows 2 and 3 are affected by MPM. After eliminating all elements in Row 1 from the“list of suspects,” Eq. 3.10 can be used to identify high-order models for ∆row2 and ∆row3 .If the polynomials B2,2 and B3,3 contain significantly non-zero entries, then ∆2,2 and ∆3,3are the MPM-affected entries.Figure 3.2: Combining S∆ and ∆ methods to pinpoint exact locations of MPM3.3 Time Delay Estimation for LTI ModelsThe identifications of S∆(q) and ∆(q) models require an adequate estimate for the timedelay. The estimation of any model requires a general input and general output, denotedas u∗(t) and y∗(t), respectively. For example, in the identification of S∆(q), u∗(t) = r(t)and y∗(t) = (t).Assume that the engineer wishes to estimate a 4th-order BJ model for both S∆(q),which would resemble the following:393.3. Time Delay Estimation for LTI ModelsAS∆(q) · y∗(t) =BS∆FS∆(q) · u∗(t) + CS∆DS∆· e(t) (3.12)Although the model is originally desired to take on a BJ structure, for the purpose ofdelay estimation it can be temporarily assumed to follow a FIR structure:y∗(t) = B(q) · u∗(t) + e(t) (3.13)where u∗(t) and y∗(t) are the general input and output time-series, respectively, foreach model. In matrix expanded form for MIMO systems, this would be expressed as:y∗1...y∗ny =B1,1 · · · B1,nu.... . ....Bny ,1 · · · Bny ,nu ·u∗1...u∗nu+e1...eny (3.14)where u∗(t) and y∗(t) represent the relevant inputs and outputs, respectively, for eachmodel. As outlined in Section 1.3, if 2N datapoints are available, then by forming theregressor:Φ(t) =u∗1(t) · · · u∗1(t−N) · · · u∗nu(t) · · · u∗nu(t−N)....... . .......u∗1(t−N) · · · u∗1(t− 2N) · · · u∗nu(t−N) · · · u∗nu(t− 2N) (3.15)and the parameter vector:403.4. Hypothesis Testing for “Significant” Model CoefficientsΘ = [B1,1,0 · · ·B1,1,N B1,2,0 · · ·B1,2,N · · · B1,nu,0 · · ·B1,nu,N ]T (3.16)as well as the time-series vectors:Y (t) = [y∗(t) · · · y∗(t−N)]T (3.17)E(t) = [e(t) · · · e(t−N)]T (3.18)the delay between input u∗(t) and output y∗(t) can be estimated by solving for Θ in:Y (t) = Φ(t) ·Θ + E(t) (3.19)This can be achieved by using any LN -norm minimization procedure for Θ [4] (the mostcommon being L1-norm and L2-norm). The number of leading near-zero terms in Θ willbe an adequate approximation of the delay.3.4 Hypothesis Testing for “Significant” Model CoefficientsDuring the second step of the MPM detection algorithm, the transfer function ∆ is tobe identified using signals u(t) and d(t). Suppose the engineer has selected a BJ modelstructure for ∆. Then the equation relating u(t) and d(t) becomes:A(q)d(t) =B(q)F (q)u(t) +C(q)D(q)e(t) (3.20)413.4. Hypothesis Testing for “Significant” Model Coefficientswhere each BJ model polynomial is expanded as:M(q) = 1 +M1 · q−1 + · · ·+MnM · q−nM (3.21)nM denotes each desired polynomial order. Assuming that the polynomials B(q) andF (q) are coprime (i.e. no common factors) and that the polynomial coefficients in F (q)are significantly larger than those in B(q), then the Bj,k(q) polynomials which contain“significantly non-zero” terms indicate contributions in matrix element ∆j,k(q) from inputsuj,k(t). Hence a corresponding MPM is present in matrix element Pˆj,k. A mathematicaldefinition of “non-zero” can be realized by making the following assumptions:1. The noise term e(t) is zero-mean, Gaussian-distributed2. As a result of 1, each polynomial coefficient Bj,k,n (n ∈ [1, · · · , nB]) is also Gaussian-distributed (but not necessarily zero-mean)Any “non-zero” polynomial coefficients in B(q) indicate the presence of MPM. However,in order to compare the polynomials against each other, a representative coefficient mustbe extracted from each polynomial. An adequate choice is the maximum absolute-valuecoefficient from each polynomial:Bi,j,max , maxn∈1,··· ,nB|Bi,j,n| (3.22)This value is treated as a mean, while its covariance is estimated using the following ap-proach outlined by [12], which involves the average squared error VN (θ, ZN ), first-derivativeof model output Ψ(t, θ), and the variance estimate λˆN :423.4. Hypothesis Testing for “Significant” Model CoefficientsVN (θ, ZN ) =1NN∑t=1122(t, θ) (3.23)Ψ(t, θ) =ddθyˆ(t|θ)θ=θˆ (3.24)λˆN =1NN∑t=12(t, θˆN ) (3.25)In this context, θ represents the polynomial coefficient estimates. Assuming that thedata size N tends to infinity, the asymptotic covariance matrix for the parameter estimatescan be evaluated as:limN→∞cov(θˆN ) = λˆN[1NN∑t=1[Ψ(t, θˆN ) ·ΨT (t, θˆN )]]−1(3.26)Finally, the covariance of each individual polynomial coefficient, cov(Bi,j,max), can besimply extracted from the diagonal entries of the matrix cov(θˆN ).To perform the hypothesis test, the first step is to construct every maximum coefficientas a normal distribution with a mean of Bi,j,max and variance cov(Bi,j,max). Next, a levelof significance α is selected by the user; standard values used in literature are 1%, 5%,or 10%. α simply represents the cumulative probability that the true value of Bi,j,maxfalls between −∞ and a threshold value on the bell curve. In other words, the true valueof Bi,j,max will fall between the threshold value and +∞ with a confidence probability of1− α.433.4. Hypothesis Testing for “Significant” Model CoefficientsFigure 3.3: Probability distribution and confidence interval of Bi,j,maxTo determine whether Bi,j,max is indeed non-zero, the user must decide whether thevalue zero lies inside or outside the confidence interval 1 − α. As indicated in Figure 3.3,only two scenarios are possible:1. Zero lies within the confidence interval 1− α, implying the probability that the truevalue Bi,j,max equals zero is “significant.”2. Zero lies outside the confidence interval 1−α, implying the probability that the truevalue Bi,j,max equals zero is “insignificant.”This basic classification can be formally constructed as a hypothesis test :• H0 (Null): Bi,j,max 6= 0, i.e. model element Pˆi,j is MPM-affected• H1 (Alternate): Bi,j,max = 0, i.e. model element Pˆi,j is MPM-freeThe test probability in this case is P (Bi,j,max = 0|∆i,j 6= 0), i.e. the probability thatthe true value Bi,j,max equals zero, given that model element Pˆi,j is MPM-affected. Therationale behind this particular definition of probability stems from the fact that if Pˆi,j443.5. Input Design for the Excitation Signalwere truly MPM-affected, then the probability that Bi,j,max = 0 should be minute, smallerthan the level of significance α. Scenario 1 corresponds to the case where P (Bi,j,max =0|∆i,j 6= 0) > α, the probability of Bi,j,max = 0 being larger than the significance level.Thus the null hypothesis is refuted and MPM is diagnosed to be absent. On the otherhand, Scenario 2 implies the case where P (Bi,j,max = 0|∆i,j 6= 0) < α, the probability ofBi,j,max = 0 being smaller than the significance level and confirming the presence of MPM.3.5 Input Design for the Excitation SignalIn typical plants, the most common types of setpoint or plant input excitations used arethe step, ramp, sinusoid, and Psuedo-Random Binary Sequence (PRBS) signals:Figure 3.4: Typical excitation signals used in industryThese signals, when injected into the system at the reference or input points, aresufficient to identify most LTI process models approximated as First-Order Plus TimeDelay (FOPTD) in the form of Pij =Kτs+1 . As stated previously in Section 3.1, theidentification of an unbiased first-order model requires at most 1 non-zero frequency in theexcitation signal. Therefore, these signals are sufficient for detecting MPM in which the453.5. Input Design for the Excitation Signalplant transfer functions (and hence the S∆(q) and ∆(q) transfer functions) are FOPTD.However, for higher-order processes, the typical excitations used for first-order processesdo not contain enough frequencies for unbiased identification. The field of Input Designjustifies the required excitation signals mathematically. In input design, a unique model isfound (for a system of any order) in a way such that the covariances of the model coefficientsare minimized. Similarly, if the excitation signal is designed to be “sufficiently exciting,”then the diagnosis of MPM using he S∆ and ∆ Sys-ID experiments outlined in Sections3.1 and 3.2 can be ascertained to be reliable. This “sufficiently exciting” signal thereforerepresents as a theoretical “lower bound” on the frequency content of the excitation signal.The starting point of input design is the definition of “persistent excitation” for SystemIdentification, which has been developed and established in [12]. The core assumptions oflinear system identification also apply here; specifically, the plant and noise models areassumed to be LTI, and the controller is assumed to be operating in a linear region. Theautocorrelation of the input signal is defined as the covariance of the input at a given time,t, with itself at a certain lag t− k, k ∈ N:Ru(k) = E[u(t) · u(t− k)] (3.27)If the input is Gaussian-distributed, it is uncorrelated with itself at all lags except zero.In this case, the autocorrelation simply equals the input covariance, since Ru(k) = 0 ∀k ∈N, except at k = 0 for which Ru(0) = E[u(t)]2. The power spectrum of the input is definedmathematically as the Discrete-Time Fourier Transform of the autocorrelation of the input:463.5. Input Design for the Excitation SignalΦu =∞∑k=−∞Ru(k) · e−jωk (3.28)Physically, the power spectrum describes how the input autocorrelation changes withany given frequency. Again, for the special case in which the input is Gaussian, Ru(0) =E[u2(t)] is the single non-zero value in the infinite sum. Therefore, the power spectrumevaluates to a constant value equal to the input covariance Ru(0) for all frequencies, re-sembling a flat line on a frequency plot:Φu =∞∑k=−∞Ru(k) · e−jωk = Ru(0) · e−jω·0 = Ru(0) = E[u(t)]2 (3.29)However, the input distribution is typically non-Gaussian, so the other terms in theinfinite sum are non-zero and the spectrum does not take on a constant value across allfrequencies.The concept of a quasi-stationary (QS) signal is required. According to [12], a signalu(t) is QS if and only if:1. E[u(t)] ≤ C,∀t ∈ N2. E[u(t)u(t− k)] ≤ C,∀k ∈ N 6= 0,∀t ∈ N3. limN→∞1N∑Nt=0 u(t) · u(t− k) ≤ C˜,∀k ∈ NC, Cˆ are arbitrary constants, and N is the data sample size. Mathematically speak-ing, the QS conditions state that the autocorrelations at different times t, as well as themean autocorrelation value averaged across all t, are bounded and finite. Now the formaldefinition of ”Persistence of Excitation” can be introduced:473.5. Input Design for the Excitation SignalTheorem 3.5.1. Persistence of Excitation (Ljung, 1999): If an input signal u(t) isQS and has a spectrum Φu(ω), then it is persistently exciting of order n if and only if, forall model transfer functions:Mn(q) = m1 · q−1 + · · ·+mn · q−nThe equality |Mn(ejw)|2 · Φu(ω) = 0 implies that Mn(ejw) = 0.The theorem states that given a particular input u(t) that is ”persistently exciting withorder n”, the nth-order model Mn(q) determined from injecting such an input is unique.The question, now, becomes how u(t) and its spectrum, Φu(ω) can be generated.Suppose that the LTI system is modelled using a BJ model structure:y(t) =B(q)F (q)· u(t) + C(q)D(q)· e(t) = G(q, θ) · u(t) +H(q, θ) · e(t) (3.30)where G(q, θ) = B(q)F (q) represents the plant model and H(q, θ) =C(q)D(q) represents thenoise model for the white noise sequence e(t). The asymptotic covariance matrix of themodel parameters, P (θN ), can be expressed as:limN→∞P (θN ) = λˆN[1NN∑t=1[Ψ(t, θˆN ) ·ΨT (t, θˆN )]]−1(3.31)where the partial derivative Ψ is taken with respect to θ for the model output yˆ:483.5. Input Design for the Excitation SignalΨ(t, θ0) ,∂∂θyˆ(t, θ0) (3.32)[10] discovered that, using Parseval’s Theorem:∫ ∞−∞|U(ω)|2dω =∞∑k=−∞|u(k)|2 (3.33)the covariance matrix P (θN ) can be re-written as a function of the input spectrum,Φu(ω), directly:P−1θ =12piλ0·∫ pi−pi[Fu(ejω, θ0) · Φu(ω) · F ∗u (ejω, θ0)]dω +12pi∫ pi−pi[Fe(ejω, θ0) · F ∗e (ejω, θ0)]dω(3.34)Conveniently, the first-order partial derivatives Fu, Fe are simply evaluated as:Fu(q, θ0) ,∂∂θ|θ=θ0G(q, θ) (3.35)Fe(q, θ0) ,∂∂θ|θ=θ0H(q, θ) (3.36)The ∗ sign denotes the complex conjugate. In order to design an adequate input forthe identification, the most straightforward choice would be to minimize the covariancematrix P (θN ), or equivalently, to maximize the information matrix P−1(θN ), and accountfor obvious physical constraints on the input spectrum. The rough optimization procedure,therefore, represents the following:493.5. Input Design for the Excitation SignalmaximizeΦu(ω)f [P−1(θN )]subject toΦu(ω) ≤ UBΦu(ω) ≤ LBwhere f represents some numerical function of the information matrix, and UB andLB some arbitrary upper and lower bounds, respectively. However, in order to numericallyevaluate the integral in 3.34, several more approximations must be made. For instance,the maximization of the information matrix can take on many forms, the simplest onesbeing the trace or eigenvalues of the matrix. Moreover, the input spectrum Φu(ω) cannotbe analyzed in its continuous form. However, it can be parameterized as a discrete sum ofbasis functions β˜k to give a discrete input spectrum with M desired coefficients [10]:Φu(ω) =M−1∑k=0ck · βk =M−1∑k=0ck · (ejωk + e−jωk) (3.37)The above equation leads to a simple numerical approximation for the function Φu(ω).Next, the physical constraints on Φu(ω) must be realized. First, Φu(ω) obviously cannotbe negative, because negative covariances do not exist. The lower bound is thereforeΦu(ω) ≥ 0. Secondly, an upper bound must exist on Φu(ω), since a rational controlengineer would not want the input covariance to be too large. This constraint is realized as12pi∫ pi−pi Φu(ω)dω ≤ α. These two constraints are difficult to formulate numerically in theirdirect forms, but an elegant approximation exists if the system falls under a special case.503.5. Input Design for the Excitation SignalSpecifically, if the input signal u(t) can be considered filtered white noise, with the filtertaking on a FIR form, then the lemma developed by [10] applies. The lemma states thatthe positivity of Φu(ω) is guaranteed if a state-space realization of the system is positive:Lemma 3.5.2. Positivity of Input Spectrum (Hjalmarsson, 2005):Let {A,B,C,D} be a controllable state-space realization of ReΦu(ω). Then, ∃Q = QTsuch that K ,Q−ATQA −ATQB−BTQA −BTQB+0 CTC D +DT ≥ 0 iff Φu(ω) ≥ 0In this special case, the K matrix can be simply evaluated as:A =O1×(M−2) 0IM−2 O(M−2)×1 , (3.38)B = [1 0 · · · 0]T , (3.39)C = [c1 · · · cM−1], (3.40)D = c0 (3.41)The upper bound of 12pi∫ pi−pi Φu(ω)dω ≤ α can be simply re-written as c0 ≤ α, wherec0 is the first coefficient in the parameterization∑M−1k=0 ck · (ejωk + e−jωk). Using thesesimplifications, the optimization procedure transforms from the search for the entire func-tion Φu(ω) into the search of coefficients ck, k ∈ 0 · · ·M − 1. One possible form of such anoptimization can be of the following, which is a convex semi-definite program (SDP) [3]:513.5. Input Design for the Excitation Signalmaxck,k∈0···M−1trace(P−1θ )subject to c0 ≤ α ⇒ 12pi∫ pi−piΦu(ω)dω ≤ αK > 0 ⇒ Φu(ω) > 0where the trace of the information matrix is maximized over all possible coefficients ck.The ck values can be transformed into the FIR filter coefficients for H(q) = h0 +h1 · q−1 +· · · + hM−1 · q−M+1 using a spectral factorization procedure developed by [16]. This willproduce a unique FIR filter H(q) (and hence a unique input u(t)) for a given spectrumΦu(ω). Finally, the actual input u(t) is then obtained by passing a white-noise sequencee˜0(t) with a specified variance through the filter H(q), i.e. u(t) = H(q) · e˜0(t).The previous section detailed the general input design procedure. An important clarifi-cation is required here: with regards to the MPM detection algorithm detailed previously,the relevant spectrum and signal are Φr(ω) and r(t), respectively. The entire input designprocedure is based on the transfer function between (t) and r(t), namely S∆(q), since the“plant” model of concern is the BJ expression in Equation 3.8. Therefore, Φr(ω) representsthe optimal setpoint spectra that minimizes the covariance of the setpoint signal, whichserves as an adequate, theoretical requirement of the frequency content in the setpoint.Formulations of the optimization problem for non-FIR cases are covered in [21]. Fur-thermore, the optimization procedure can be refined further by the addition of Markovchains [5]. Using this approach, the engineer can limit the frequency band of the optimalinput signal within a desired range, if the frequencies of the optimal input generated fromthe previously mentioned method are too high. Since the optimization problem is altered,the optimal input signal will be different, and the optimization problem itself may have no523.6. Computational Savings of the Proposed Methodfeasible solutions. Nevertheless, it is a necessary refinement for plants that cannot sufferhigh-frequency excitation.3.6 Computational Savings of the Proposed MethodConsider the general MIMO plant with nu inputs and ny outputs, as depicted below:P =P11 · · · P1nu.... . ....Pny1 · · · Pnynu (3.42)Assume the ”worst-case” scenario where MPM exists in every row. The table belowcompares the number of setpoint excitations required for identification, as well as thenumber of models to be identified, for the traditional case of plant re-identification and theproposed method, respectively:Table 3.1: Computational Requirement Comparison# of excitations # of polynomials to be IDedTraditional plant re-ID nu ny · nuProposed algorithm ny n2y + nR · nuThe excitation column implies that the wider the plant matrix (ny < nu), the lessnumber of setpoint excitations required for MPM detection, as compared to traditionalplant re-identification. This has the direct consequence of reduced interference on normalplant operation, since even simple bump tests can be extremely disruptive. On the otherhand, square plants (ny = nu) offer no reduction in required excitations, and tall plants(ny > nu) require more excitation than traditional re-identification.533.6. Computational Savings of the Proposed MethodIn regards to the identification efforts involved, again consider the case where the plantmatrix is wide (ny < nu). In this case, n2y < ny · nu, and if MPM is sparse, then nR ≈ 0,and therefore n2y+nR ·nu < ny ·nu, resulting in a significant reduction in identified models.On the other hand, for square matrices sparsely affected by MPM, ny = nu, n2y = ny · nu,and n2y +nR ·nu > ny ·nu, resulting in a slight increase in identified models. Finally, if theplant matrix is tall (ny > nu), then n2y > ny · nu and therefore n2y + nR · nu >>> ny · nu,resulting in a significant increase in identified models. In light of these results, it mustbe stressed again that the proposed algorithm only provides significant reductions in bothexcitations and computational savings if MPM is sparse and the plant matrix is wide.Example 3.6.1. This simple example demonstrates the reductions achieved in both exci-tations and models to be identified, using the proposed algorithm.Consider a plant with 100 inputs and 50 outputs (ny = 50 and nu = 100), sufferingfrom sparse MPM (only 5 of its rows are MPM-affected, i.e. nR = 5). Using Table ??,the number of required excitations and identified models are as follows:Table 3.2: Computational Requirement Comparison# of excitations # of polynomials to be IDedTraditional plant re-ID 100 5000Proposed algorithm 50 3000The proposed method requires only half the number of excitations and slightly morethan half the number of identified models, compared to traditional plant re-identification.Despite the possible computational savings achievable, specific re-tuning of controllerelements is only recommended if the proposed algorithm detects sparse MPM. As demon-543.6. Computational Savings of the Proposed Methodstrated numerous times previously, if the plant is severely affected by MPM (with nR ≈ ny),then the control engineer might as well perform plant re-identification rather than use theproposed MPM detection algorithm.55Chapter 4Simulation ResultsIn order to assess the feasibility and robustness of the previously proposed algorithm forMPM detection, the algorithm is submitted to a series of artificial exercises. First, the S∆test is conducted on a 3x3 plant containing no MPM, to demonstrate the absence of anyfalse diagnosis of MPM. Then, the combined S∆ and ∆ test is conducted on a 3x3 plantcontaining sparse MPM entries, to demonstrate its ability to correctly pinpoint the MPM-affected elements. The last simulation involves the MPM detection of a 5x5 plant, whichalso contains sparse MPM entries. This shows that the algorithm retains its robustnesswhen tested on larger systems. Finally, the algorithm is applied to closed-loop plant datafrom an industrial lime slaker. Altogether, the algorithm is shown to be reliable whentested on any system that follows the assumptions outlined in Section 3.4.1 3x3 Plant: Case of No MPMThe feasibility of the proposed algorithm is first demonstrated on a 3x3 (nu = ny = 3)plant artificially constructed in Matlab. The plant uses similar transfer functions as thatof the famous Shell Benchmark problem, with reduced delays:564.1. 3x3 Plant: Case of No MPMP =4.05·e−6s50s+11.77·e−7s60s+15.88·e−6s50s+15.39·e−4s50s+15.72·e−3s60s+16.9·e−3s40s+14.38·e−5s33s+14.42·e−5s44s+17.219s+1 (4.1)In this case, no MPM is present, therefore Pˆ = P . The goal is to observe whether thealgorithm falsely identifies the presence of MPM, when none actually exists.The MPC controller has the following tuning parameters:• Sampling time: 1 sec• Prediction Horizon: 100 sec• Control Horizon: 30 sec• Input Weights: diag [1 1 1]• Input Rate Weights: diag [0.1 0.1 0.1]• Output Weights: diag [1 1 1]• Input Constraints: None• Input Rate Constraints: None• Output Constraints: NoneThe rest of the relevant control signals are defined as follows:• Output: corrupted with two types of disturbances:1. Measurement Noise: zero-mean white noise with Standard Deviation (SD) =0.005574.2. 3x3 Plant: Case of Sparse MPM2. Unmeasured Disturbance: zero-mean white noise with SD = 0.0075, throughfirst-order filter diag [ zz−0.95zz−0.95zz−0.95 ]• Setpoint (Reference): Psuedo-Random Binary Sequence (PRBS) with frequencyband [0 fN ] (fN , Nyquist frequency)First, the S∆ matrix is identified using Matlab’s “bj” function, with each row of thematrix having 4th-order BJ polynomials (i.e. nA = nB = nC = nD = 4). The B polynomi-als contain delays of [nk,B1 nk,B2 nk,B3 ], estimated by Matlab’s “delayest” function. Then,the singular values of S∆ are plotted for each row in the S∆ matrix. Observe that since noMPM is present, both the overall and row ||S∆(ω)||∞ values are reasonably close to zero.Figure 4.1: Overall S∆(ω) plot for the no-MPM caseSince the S∆ test shows that no MPM is present anywhere within the plant, ∆ identi-fication need not be performed.4.2 3x3 Plant: Case of Sparse MPMThe same MPC controller and tuning parameters (as described in the previous section) areused. In this case, however, modelling errors are artificially introduced on transfer matrixelements Pˆ12, Pˆ13, and Pˆ21. Specifically, while the plant transfer matrix remains as:584.2. 3x3 Plant: Case of Sparse MPMP =4.05·e−6s50s+11.77·e−7s60s+15.88·e−6s50s+15.39·e−4s50s+15.72·e−3s60s+16.9·e−3s40s+14.38·e−5s33s+14.42·e−5s44s+17.219s+1 (4.2)The plant model matrix now becomes affected by sparse MPM (affected elements indi-cated in red):Pˆ =4.05·e−6s50s+11.77·e−4s45s+13.5·e−6s50s+15.39·e−4s70s+15.72·e−3s60s+16.9·e−3s40s+14.38·e−5s33s+14.42·e−5s44s+17.219s+1 (4.3)In the MPM (∆) matrix, ∆12 is a combination of time constant and delay mismatch,∆13 a severe gain mismatch, and ∆21 a time constant mismatch. After applying the S∆test, non-unity values of ||S∆(ω)||∞ are observed in Rows 1 and 2, indicating the presenceof MPM for these rows. For Row 3, ||S∆(ω)||∞ = 1, meaning that this row is MPM-free.Figure 4.2: Overall and row S∆(ω) plots for the sparse MPM case594.2. 3x3 Plant: Case of Sparse MPMThe question now is whether the MPM elements ∆12,∆13, and ∆21 can be correctlypinpointed. To confirm this, the ∆ test is conducted by identifying models between d1and [u1 u2 u3] and between d2 and [u1 u2 u3], in the form of 4th order BJ polynomials(nA = nB = nC = nD4). The coefficients for the two relevant polynomials, B and F , areshown below. A delay of k units is represented by a multiplication of each coefficient withq−k (q , forward-shift operator):Table 4.1: BJ model between d1 and uPolynomial q−1 q−2 q−3 q−4 Delay Fit(%)Bu1 −2.084 · 10−4 6.264 · 10−4 1.047 · 10−4 4.149 · 10−4 7 99.37Bu2 −3.923 · 10−2 −1.595 · 10−3 −8.174 · 10−4 2.856 · 10−2 5Bu3 −6.525 · 10−4 −3.240 · 10−4 4.667 · 10−2 9.816 · 10−4 5Fu1 1.000 −1.833 9.761 · 10−2 1.433Fu2 1.000 −2.345 1.857 5.254 · 10−1Fu3 1.000 2.051 · 10−1 −1.450 −3.388 · 10−1Table 4.2: BJ model between d2 and uPolynomial q−1 q−2 q−3 q−4 Delay Fit(%)Bu1 −6.943 · 10−3 −1.945 · 10−4 1.300 · 10−4 −5.290 · 10−4 5 99.68Bu2 6.211 · 10−4 2.838 · 10−4 2.922 · 10−4 7.952 · 10−4 5Bu3 6.971 · 10−4 1.867 · 10−4 −1.914 · 10−4 2.671 · 10−4 2Fu1 1.000 −1.464 8.842 · 10−1 7.499 · 10−1Fu2 1.000 −2.938 3.285 −1.694Fu3 1.000 −1.689 3.213 · 10−1 7.900 · 10−1In each table, “significantly non-zero” coefficients have been highlighted in red. The BJmodel fits to estimation data are reasonably accurate (> 99%). In the BJ model between d1and u (Table 4.1), the polynomials Bu2 and Bu3 contain significantly non-zero coefficients,compared to all other coefficients. Similarly, in the model between d2 and u (Table 4.2),the polynomial Bu1 contains a significantly non-zero coefficient compared to all others.These results indicate that CV1 is suffering from MPM due to inputs u2 and u3, and CV2is suffering from MPM due to input u1.604.2. 3x3 Plant: Case of Sparse MPMThe following table consolidates the Bi,j,max coefficients for each BJ polynomial. Foreach maximum coefficient, a hypothesis test was conducted to determine whether it couldbe considered “non-zero.” For coefficients B1,2,max, B1,3,max, and B2,1,max, the probabilityof their true values equating to zero were found to be lower than the user-defined 1%significance level, thus confirming the presence of MPM in model elements Pˆ12, Pˆ13, andPˆ21. On the other hand, for B1,1,max, B2,2,max, and B2,3,max, the probability of equatingto zero were found to be higher than the significance level of 1%. This indicates that thesecoefficients cannot be distinguished from zero, and hence the corresponding model elementsPˆ1,1, Pˆ2,2, and Pˆ2,3 are likely to be MPM-free.Table 4.3: Bi,j,max coefficients for 3x3 sparse MPM caseBi,1,max Bi,2,max Bi,3,maxd1 6.264 · 10−4 3.923 · 10−2 4.667 · 10−2d2 6.943 · 10−3 7.952 · 10−4 6.971 · 10−4d3Table 4.4: Hypothesis testing for coefficients for 3x3 sparse MPM case (α = 1%)P (Bi,1,max = 0)(%) P (Bi,2,max = 0)(%) P (Bi,3,max = 0)(%)d1 2.85 0.00 0.00d2 0.13 5.90 6.62d3In order to refine the excitation signals, the input design procedure from Section 3.5 isused. For example, using the available S∆ data, the optimal setpoint spectra and excitationsignal can be determined. Using an optimal setpoint length of 100 sampling intervals, anupper constraint of α = 10 on the setpoint spectrum integral, and a noise variance of0.0075, the optimal results are:614.3. 5x5 Plant: Case of Sparse MPMFigure 4.3: Optimal spectrum and setpoint signals for MPM DetectionThe optimal spectrum for Setpoint 1 contains a high-frequency component at 3 rads ,and the optimal excitation signal resembles a PRBS signal with varying amplitude. On theother hand, the optimal spectrum for Setpoint 2 contains mostly low-frequency componentsbelow 10−1 rads , and hence the corresponding optimal excitation resembles a low-frequencysinusoid. Since the exact noise variance is known from simulation, the optimization proce-dure gives (one of many) correct possible solutions for the optimal setpoint.4.3 5x5 Plant: Case of Sparse MPMThe final exercise involves a 5x5 plant which suffers from sparse MPM. Specifically, theplant suffers from modelling errors in elements Pˆ11, Pˆ21, Pˆ23, Pˆ42, and Pˆ45. Suppose all5x5=25 elements of the true plant transfer matrix are 5·e−5s50s+1 , without loss of generality. Themodel matrix contains following elements (with the MPM-affected elements highlighted inred):624.3. 5x5 Plant: Case of Sparse MPMPˆ =5·e−8s50s+15·e−5s50s+15·e−5s50s+15·e−5s50s+15·e−5s50s+14·e−5s60s+15·e−5s50s+15·e−5s25s+15·e−5s50s+15·e−5s50s+15·e−5s50s+15·e−5s50s+15·e−5s50s+15·e−5s50s+15·e−5s50s+15·e−5s50s+13.5·e−5s50s+15·e−5s50s+15·e−5s50s+15·e−3s75s+15·e−5s50s+15·e−5s50s+15·e−5s50s+15·e−5s50s+15·e−5s50s+1(4.4)Similar MPC controllers and reference signals are constructed as in the previous sub-section. To detect MPM, the S∆ test is first performed on the plant. Non-unity valuesof ||S∆(ω)||∞ are observed in Rows 1, 2, and 4, which indicate the correct diagnosis ofMPM in these rows. In Rows 1 and 2, the MPMs consist of severe time constant and delaymismatches, while the gain mismatches are mild. Therefore the peaks of S∆(ω) are smallerin these rows, compared to Row 4 where a severe gain mismatch exists in element ∆42.Rows 3 and 5 show ||S∆(ω)||∞ values of 1, correctly diagnosing the absence of MPM inthese rows.634.3. 5x5 Plant: Case of Sparse MPMFigure 4.4: Overall and row S∆(ω) plots for a 5x5 plant suffering from sparse MPM.Next, ∆ is identified between each MPM-affected simulation errors d1, d2, and d4, andinputs [u1 · · ·u5]. The 4th-order BJ coefficients for polynomials B and F are tabulatedbelow, again with the significantly non-zero elements highlighted in red:644.3. 5x5 Plant: Case of Sparse MPMTable 4.5: 5x5 plant: BJ model between d1 and uPolynomial q−1 q−2 q−3 q−4 Delay Fit(%)Bu1 1.003 9.623 · 10−1 9.201 · 10−3 −8.009 · 10−1 6 99.98Bu2 2.301 · 10−2 3.580 · 10−1 2.357 · 10−1 2.934 · 10−2 0Bu3 −1.100 · 10−3 2.444 · 10−4 −1.390 · 10−2 −3.256 · 10−3 0Bu4 −1.340 · 10−3 −2.464 · 10−3 3.918 · 10−3 −3.396 · 10−2 0Bu5 −2.638 · 10−3 −1.950 · 10−2 2.637 · 10−2 −2.606 · 10−2 0Fu1 1.000 9.907 · 10−1 9.232 · 10−3 −8.104 · 10−1Fu2 1.000 4.968 · 10−1 −2.260 · 10−2 −6.002 · 10−1Fu3 1.000 −5.471 · 10−1 −3.288 · 10−1 −2.937 · 10−1Fu4 1.000 −1.776 · 10−1 −1.023 3.177 · 10−1Fu5 1.000 −2.263 2.693 −1.952Table 4.6: 5x5 plant: BJ model between d2 and uPolynomial q−1 q−2 q−3 q−4 Delay Fit(%)Bu1 1.103 · 10−2 −3.501 · 10−3 −4.371 · 10−2 −9.448 · 10−2 5 97.52Bu2 −2.338 · 10−3 3.585 · 10−2 2.351 · 10−2 2.934 · 10−2 0Bu3 −1.111 · 10−3 2.470 · 10−4 −1.391 · 10−1 −3.322 · 10−2 5Bu4 −1.307 · 10−3 −2.549 · 10−3 3.909 · 10−3 3.961 · 10−3 0Bu5 −2.600 · 10−3 −1.597 · 10−2 2.636 · 10−2 −2.608 · 10−2 0Fu1 1.000 −7.624 · 10−1 5.136 · 10−1 −9.998 · 10−1Fu2 1.000 1.404 7.790 · 10−2 7.207 · 10−1Fu3 1.000 −2.680 3.107 −1.838 · 10−1Fu4 1.000 −5.201 · 10−1 −3.535 · 10−1 −1.633 · 10−1Fu5 1.000 −3.181 3.746 −1.942654.3. 5x5 Plant: Case of Sparse MPMTable 4.7: 5x5 plant: BJ model between d4 and uPolynomial q−1 q−2 q−3 q−4 Delay Fit(%)Bu1 −1.510 · 10−2 −2.873 · 10−2 −1.694 · 10−2 5.797 · 10−2 0 99.62Bu2 9.003 · 10−3 −1.677 · 10−2 9.778 · 10−3 1.339 · 10−1 4Bu3 1.149 · 10−3 5.662 · 10−4 7.836 · 10−2 −3.555 · 10−3 0Bu4 −1.480 · 10−2 −2.935 · 10−3 4.317 · 10−2 −6.888 · 10−3 0Bu5 9.991 · 10−1 −9.737 · 10−1 −5.231 · 10−1 2.473 · 10−1 4Fu1 1.000 −2.633 2.302 −6.780 · 10−1Fu2 1.000 −2.434 2.021 · 10−1 −5.416 · 10−1Fu3 1.000 −1.425 9.320 · 10−1 −3.452 · 10−1Fu4 1.000 −5.276 · 10−1 −9.489 · 10−1 5.523 · 10−1Fu5 1.000 −9.736 · 10−1 −6.237 · 10−1 2.422 · 10−1In the BJ model between d1 and [u1 · · ·u5] (Table 4.5), the polynomials Bu1 containsignificantly non-zero coefficients, compared to all others. This is also the case in the modelbetween d2 and [u1 · · ·u5] (Table 4.6) for the polynomials Bu1 and Bu3 , and in the modelbetween d4 and [u1 · · ·u5] (Table 4.7) for the polynomials Bu2 and Bu5 . The tables belowshow that for the max coefficients B1,1,max, B2,1,max, B2,3,max, B4,2,max, and B4,5,max, theprobability of equating to zero is less than the significance level of 1%, thus confirmingthat the model elements Pˆ1,1, Pˆ2,1, Pˆ2,3, Pˆ4,2, and Pˆ4,5 are indeed affected by MPM. Onthe other hand, max coefficients B1,2,max, B1,3,max, B1,4,max, B1,5,max, B2,2,max, B2,4,max,B2,5,max, B4,1,max, B4,3,max, and B4,4,max have a likelihood of equating to zero greater than1%, therefore indicating that these elements are likely to be MPM-free. Therefore, thehypothesis test has correctly identified the MPM-affected elements, namely ∆11, ∆21, ∆23,∆42, and ∆45 in the original, 5x5 artificial plant.664.3. 5x5 Plant: Case of Sparse MPMTable 4.8: Bi,j,max coefficients for 5x5 sparse MPM case|Bu1|max |Bu2|max |Bu3|max |Bu4|max |Bu5|maxd1 1.003 3.580 · 10−1 1.390 · 10−2 3.396 · 10−2 2.637 · 10−2d2 9.448 · 10−2 3.585 · 10−2 1.391 · 10−1 3.961 · 10−3 2.636 · 10−2d3d4 5.797 · 10−2 1.339 · 10−1 7.836 · 10−2 4.317 · 10−2 9.991 · 10−1d5Table 4.9: Hypothesis testing for coefficients for 5x5 sparse MPM case (α = 1%)P (|Bu1|max = 0)(%) P (|Bu2|max = 0)(%) P (|Bu3|max = 0)(%) P (|Bu4|max = 0)(%) P (|Bu5|max = 0)(%)d1 0.00 4.73 5.12 5.64 3.38d2 0.00 10.01 0.00 3.15 3.76d3d4 4.99 0.00 5.01 4.84 0.00d5Using the input design procedure from Section 3.5 and the available S∆ data, theoptimal setpoint spectra and excitation signal are determined. The parameters for theoptimization procedure are an optimal setpoint length of 100 sampling intervals, an upperconstraint of α = 10 on the setpoint spectrum integral, and a noise variance of 0.0075:674.3. 5x5 Plant: Case of Sparse MPMFigure 4.5: Optimal spectrum and setpoint signals for MPM DetectionAll optimal spectra contain a high-frequency components at or above 1 rads , and theoptimal excitation signals resemble PRBS signals with varying amplitudes.68Chapter 5Application of MPM DetectionAlgorithm on Industrial PlantDataIn order to demonstrate the effectiveness of the proposed MPM detection algorithm on real,closed-loop plant data collected from industry, the particular example of the papermakingprocess is used. The most popular form of chemical pulping is known as the “Kraft process,”which converts wood chips into usable wood pulp by dissolving the chemical bonds thatlink the lignin to cellulose. By doing so, the wood pulp consists mostly of cellulose fibers,which allows it to be easily treated and pressed into paper sheets of in the further stagesof the process.The typical Kraft process is shown in the following diagram. The raw wood chipsare first pre-treated with steam and filled with white liquor (an alkali solution consistingmainly ofNa2(OH) andNa2S), in order to remove the air within and fill the chip capillarieswith the appropriate chemicals for ”cooking”. In the ”cooking” step, the chips are tossedinto vessels known as ”digesters” at 180 − 190◦C and 2 − 3 atm, where the white liquorslowly dissolves the chemical bonds between the lignin/ hemicellulose particles and thecellulose fiber structures. As the cooking step progresses, a brownish-black liquid knownas black liquor is formed, which is a mixture of lignin residue, hemicellulose waste, and69Chapter 5. Application of MPM Detection Algorithm on Industrial Plant Dataother inorganic materials removed from the wood chips.Figure 5.1: Kraft process for the delignification of wood chipsAfter the chips are cooked, on one hand, the resulting wood pulp is collected andwashed to remove any lingering chemicals. On the other hand, the black liquor is collectedand concentrated in a multiple-stage evaporator, turning the liquor from a wet 14-18 %solids composition into a drier 65-75 % solids composition. Then, the resulting black liquoris combusted in a boiler in order to recover the necessary NaOH, Na2S, and Ca(OH)2used for the delignification process. This process is known as “re-causticizing.” First, theleftover Na2SO4 is reduced to Na2S using organic material available in the black liquor:Na2SO4 + 2C −→ Na2S + 2CO2 (5.1)Next, green liquor, a mixture of Na2S and Na2CO3, is reacted with Ca(OH)2 in the70Chapter 5. Application of MPM Detection Algorithm on Industrial Plant Datacausticizing reaction:Na2CO3 + Ca(OH)2 ←→ Na2S + 2NaOH + CaCO3 (5.2)In order to regenerate the Ca(OH)2 used in the previous reaction, the resulting CaCO3is heated in a lime kiln:CaCO3 −→ CaO + CO2 (5.3)Finally, the CaO is reacted with water to produce the required Ca(OH)2 in the slakingreaction:CaO +H2O −→ Ca(OH)2 (5.4)Figure 5.2: Lime slaker for the conversion of lime to calcium hydroxide.71Chapter 5. Application of MPM Detection Algorithm on Industrial Plant DataThe practical control problem at hand revolves around this last chemical reaction,which occurs in a lime slaker, a large vessel with a hopper at one end for the insertionof the raw lime CaO, and a chute at the other end for the output of Ca(OH)2. Theproduced Ca(OH)2 is subsequently consumed in the causticizing reaction to form CaCO3,which occurs in parallel with the slaking reaction [19]. The middle section of the slakerconsists of a reactor, which thoroughly stirs the calcium oxide and water, and containsvarious temperature and pressure controllers. Temperature is the key control variable; ifthe temperature is too high, the water within the slaker mixture will boil, causing twosafety hazards. First, the boiling water may spill out of the vessels and onto the floorof the plant. Secondly, a high temperature encourages the precipitation of lime (CaO),especially if excess lime is present in the inlet. Precipitation causes clogging of white liquorfilters and loss of Ca(OH)2 product. The following temperature-related control variables,strategies, and their rationales, are as follows:• Temperature Rise (y1): the temperature difference between the slaker inlet and outlet,which has a hard physical constraint. Controlled by manipulating the lime-to-liquorinlet flow ratio. Decreasing the amount of lime present lessens the temperature rise,since the total enthalpy of reaction is decreased, and vice versa. Decreasing the slakerinlet temperature while holding the lime-to-liquor ratio constant, however, does notaffect the temperature rise, since the total enthalpy of reaction remains the same.• Slaker outlet temperature (y2): manipulated by both the inlet flow temperature andlime-to-liquor flow ratio. Decreasing the inlet flow temperature with the same lime-to-liquor ratio decreases the product temperature, due to the same enthalpy of reaction.Decreasing the lime-to-liquor temperature with the same inlet temperature decreasesthe enthalpy of reaction, thus also decreasing the product temperature.72Chapter 5. Application of MPM Detection Algorithm on Industrial Plant DataAs mentioned previously, the main manipulated variables are the lime-to-liquor ratio(u1), and inlet feed temperature (u2). The main control variables are the slaker temper-ature rise (y1), and the slaker outlet temperature (y2). After normalizing all variables indeviation form, the dynamic relationship between these variables can be approximated asthe following FOPTD Laplace-domain transfer functions:P = 1700s+1 · e−120s 00.6700s+1 · e−120s 2300s+1 · e−90s (5.5)18 hours worth of operational plant data are available for these four variables, eachsampled at intervals of 5 seconds:Figure 5.3: Normalized process variables for the lime slakerSetpoint excitation, r1, is available for the normalized slaker temperature rise y1 in theform of step changes. r1 is plotted with y1 to show the system’s dynamic response to the73Chapter 5. Application of MPM Detection Algorithm on Industrial Plant Datasetpoint changes:Figure 5.4: Normalized slaker temperature rise y1 (red) and slaker temperature rise setpointr1 (black)y2, the normalized slaker outlet temperature, does not have a setpoint; it simply hasa high constraint of 0.62. The strict constraint exists to prevent the boiling of the slakermixture and precipitation of lime, as mentioned previously. Therefore, r2 does not existand no setpoint excitation is available for CV2. Despite the lack of setpoint excitation inboth CVs (compared to a theoretical white-noise or PRBS signal used for excitation), theresults from the S∆ test are as follows:74Chapter 5. Application of MPM Detection Algorithm on Industrial Plant DataFigure 5.5: ||S∆||∞ values for the lime slakerSince ||S∆||∞ > 1 for both CV1 and CV2, MPM is suspected to exist in both variables.In order to pinpoint the exact locations of MPM, the MPM matrix ∆ is identified using4th-order BJ models. The maximum absolute-values of the B polynomials are as follows;element 12 is omitted, due to the physical considerations explained previously, such thatP12 = Pˆ12 = 0:Table 5.1: Bi,j,max coefficients for industrial lime slakerBi,1,max Bi,2,maxd1 8.10 · 10−3 -d2 1.21 · 10−2 6.51 · 10−2The standard deviations of each Bi,j,max coefficient are:75Chapter 5. Application of MPM Detection Algorithm on Industrial Plant DataTable 5.2: Standard deviations for Bi,j,max coefficients for industrial lime slakerBi,1,max Bi,2,maxd1 6.96 · 10−4 -d2 2.03 · 10−5 4.68 · 10−5Finally, using the assuming that each Bi,j,max coefficient is Gaussian-distributed, theprobabilities that each Bi,j,max equals zero are:Table 5.3: Hypothesis testing for coefficients for industrial lime slakerP (Bi,1,max = 0)(%) P (Bi,1,max = 0)(%)d1 0.00 -d2 0.00 0.00Every probability is close to zero, therefore even if a strict test statistic of α = 1% isused, the hypothesis test firmly concludes that B1,1,max 6= 0,B2,1,max 6= 0, and B2,2,max 6= 0,and hence MPM exists in the corresponding elements in the model matrix.If the control engineer wished to refine the excitation signals used such that the S∆ and∆ models used for MPM detection could be firmly ascertained, an iterative input designapproach could be employed. The data available from the first pass, namely the S∆ modelcoefficients, can be subject to the input design experiment detailed in Section 3.5 in orderto find the optimal setpoint spectra Φr1(ω), Φr2(ω), and hence optimal spectrum signalsr1(t), r2(t). An example set of the optimal setpoint spectra and signals, with a length of100 sampling intervals, an assumed noise variance of 0.01, and α = 10 is:76Chapter 5. Application of MPM Detection Algorithm on Industrial Plant DataFigure 5.6: Optimal spectrum and setpoint signals for MPM DetectionNote that both optimal setpoint signals contain high-frequency components, and re-semble a sum of PRBS signals. The control engineer must then, according to the processat hand, decide whether these optimal excitation signals are acceptable for the process. Ifthey are acceptable, then the control engineer can excite the process using these signalsand repeat the MPM detection algorithm. If they are unacceptable, then the engineer cansimply stop here and accept the quality of the MPM detection achieved thus far.The (assumed) noise variance heavily affects the optimal excitation. In this example,an arbitrary but reasonable noise variance was assumed, since the process contained un-measured disturbances. In most industrial plants approximated by LTI and/or FOPTDsystems, the disturbance can be assumed to be white or first-order filtered, as long as adrifting disturbance does not exist. The characteristics of the disturbance and its corre-sponding variance must be approximated by making educated guesses from physically orcontrol-based intuitions, if they cannot be discerned directly from closed-loop plant data.77Chapter 6Conclusions6.1 Summary of Present ContributionsThis thesis has detailed the mathematical rationale behind the proposed 2-step algorithm.The first step detects the existence of MPM by identifying the S∆ matrix between thecontroller inputs and setpoints r; the maximum singular value of this matrix is equivalentto the H∞ norm of the ratio between the actual and designed control errors. The secondstep pinpoints the exact elements in the model matrix that are affected by MPM. Inthe second step, the MPM-affected rows of the plant matrix are singled out, and theadditive mismatch matrix ∆ is identified. Specifically, high-order models are identifiedbetween simulation errors d and all inputs u. The inputs which correspond to the non-zeropolynomial coefficients in the (∆) matrix are the specific ones causing MPM. The “non-zero” coefficients are ascertained using hypothesis testing, and optimal setpoint excitationsignals are computed using input design should the engineer feel the need to repeat theentire identification procedure.The most important underlying assumptions of the proposed algorithm are that thesystem is approximately LTI within the relevant range of operation, the setpoint is “per-sistently exciting,” and that the disturbance is zero or first-order filtered white noise.Moreover, the MPM must be sparsely-distributed throughout the plant in order for thealgorithm to retain its computational superiority compared to traditional re-identification.786.2. Areas for Future WorkThe effectiveness of this algorithm was demonstrated on artificial 3x3 and 5x5 plants suf-fering from sparse MPM, as well as a 2x2 industrial lime slaker model. In every case, theS∆ test has correctly identified the MPM-affected CVs. The ∆ identification step, aidedby hypothesis testing, has correctly pinpointed the exact MPM-affected elements in everyMPM-affected row of the plant matrix. Finally, the optimal setpoint excitation sequencefor every case has been determined using the available closed-loop data and methods frominput design.6.2 Areas for Future WorkIn this final section, the main theoretical and practical limitations of the proposed MPMdetection algorithm are explored. These limitations are considered first from the perspec-tive of Sys-ID alone, then from the more practical perspective of control engineers and/orcontrol software, which may possess limited control-related analytic capabilities.The first two important assumptions are sufficient setpoint excitation and zero or first-order filtered noise. If any of these assumptions are not fulfilled, the proposed algorithmmay fail to detect MPM reliably. The results from Sections 4.2, 4.3, and 5 have shown thatthe optimal setpoint excitations are high-frequency signals, which may be unacceptable inpractice. Specifically, if these optimal excitations are not used, false positive and/or missednegatives of MPM detections may be observed, since the S∆(q) and ∆(q) matrices may be-come biased and/or non-unique. [8] and [1] have both suggested MPM detection approachesusing partial correlations between simulation errors d and inputs u. For these methods,the amount of excitation required does not exceed that provided by typical bump tests.Therefore, these approaches can be combined with the proposed Sys-ID-based method inorder to minimize the amount of excitation required. Moreover, the presence of non-zeromodel coefficients in the ∆(q) matrix are used to detect whether the simulation errors d(t)796.2. Areas for Future Workcontains a contribution from an input uj(t), j ∈ 1, · · · , nu. Currently, the hypothesis testmethod used in Section 3.4 assume Gaussian-distributed model parameters. If the noiseis no longer white or first-order, then the hypothesis test may yield inconsistent results.The same can be said for the input design procedure, where if the nature of disturbanceis unknown, the noise variance used to calculate the optimal spectrum may be inaccurate,and the calculated optimal will be incorrect compared to the true optimal. Moreover, theconclusions of the hypothesis test depend on the threshold probability α chosen by the en-gineer. The selection of α requires a more mathematically-sound justification, but this maybe difficult to design. A viable alternative could possibly be achieved by using conceptsin alarm management, where false-detection and missed-detection rates [9] are providedalong with the conclusions of the hypothesis test with a given α.Practically speaking, the proposed method provides computational savings compared toplant re-identification, only when the plant matrix is wide (i.e. more inputs than outputs).The computational complexities of the hypothesis test and input design are not consid-ered, since they vary with the analytic software available to the plant. The probabilisticand frequency-domain calculations may be too burdensome to handle in conventional con-trol software. In these cases, the control engineer must be satisfied with the quality ofsetpoint excitation available (usually in the form of step changes), because the optimalhigh-frequency excitation signals are unknown and/or cannot be used. To improve theinformation provided by the MPM algorithm, the type of MPM (gain, time constant, timedelay), as well as the relative severity of each type in comparison to one another, shouldalso be made available at the end of the diagnosis. If these are known, the control engineercould save further time and re-identification effort. For instance, delay and time constantmismatches are usually more malignant compared to gain mismatches; model gains canusually be changed incrementally to recover controller performance, but delay and time806.2. Areas for Future Workconstants cannot. Additionally, the type and severity indices also provide physical insightto the causes behind the observed MPM. For example, gain MPMs are usually associatedwith process parameter changes over time, while the process remains in a linear range ofoperation [15]. On the other hand, time constant and delay MPMs indicate process dy-namic changes over time, which are significantly more difficult to rectify. A severity indexcomparison would allow an educated decision on whether the re-tuning of few transferelements is feasible, as opposed to re-identifying the entire plant model.81Bibliography[1] A. Badwe, R. Gudi, and R. Patwardhan. Detection of model-plant mismatch in mpcapplications. Journal of Process Control, 19:1305–1313, 2009.[2] A. Badwe, R. Patwardhan, and S. Shah. Quantifying the impact of model-plantmismatch on controller performance. Journal of Process Control, 20:408–425, 2010.[3] S. Boyd and L. Vandenberghe. Semidefinite programming relaxations of non-convexproblems in control and combinatorial optimization. Communications, Computation,Control and Signal Processing, 15:279–288, 1997.[4] S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press,Cambridge, 2004.[5] C. Brighenti, B. Wahlberg, and C. R. Rojas. Input design using markov chains forsystem identification. Joint 48th IEEE Conference on Decision and Control and 28thChinese Control Conference, pages 1557–1562, 2009.[6] R. Carlsson. A practical approach to detection of plant model mismatch for mpc.2010.[7] H. S. Fogler. Elements of Chemical Reaction Engineering, Fourth Edition. Pearson,New Jersey, 2006.82Bibliography[8] N. Iqbal, N. Yusoff, and L. Tufa. Comparison between arx and fir decorrelationmethods in detecting model-plant mismatch. Journal of Applied Sciences 14, 15:1711–1719, 2014.[9] I. Izadi, S. Shah, D. Shook, S. Kondaveeti, and T. Chen. Optimal alarm design.Proceedings of 7th IFAC Symposium on fault detection, supervision and safety oftechnical processes, pages 645–650, 2009.[10] H. Jansson and H. Hjalmarsson. Input design via lmis admitting frequency-wisemodel specifications in confidence regions. IEEE Transactions on Automatic Control,50:1534–1549, 2005.[11] M. Kano, Y. Shigi, and S. Hasabe. Detection of significant model-plant mismatch fromroutine operation data of model predictive control system. International Symposiumon Dynamics and Control of Process Systems, 2010.[12] L. Ljung. System Identification: Theory for the User. Prentice Hall, New Jersey, 1999.[13] J. Maciejowski. Predictive Control with Constraints. Prentice Hall, London, 2002.[14] K. S. McClure. Algorithms for nonlinear process monitoring and controller perfor-mance recovery with an application to semi-autegenous grinding. 2013.[15] L. Olivier and I. Craig. Model-plant mismatch detection and model update for a run-of-mine ore milling circuit under model predictive control. Journal of Process Control,23:100–107, 2013.[16] A. Papoulis. Signal Analysis. McGraw-Hill, Michigan, 1977.[17] D. E. Seborg, T. F. Edgar, D. A. Mellichamp, and F. J. Doyle. Process Dynamics andControl, Third Edition. Wiley, New Jersey, 2011.83Bibliography[18] S. Skogestad and I. Postlethwaite. Multivariable Feedback Control, Second Edition.Wiley, London, 2005.[19] H. Tran and E. K. Vakkilainnen. The kraft chemical recovery process. IncreasingEnergy and Chemical Recovery Efficiency in the Kraft Process, pages 1–8, 2008.[20] J. Webber and Y. Gupta. A closed-loop cross-correlation method for detecting modelmismatch in mimo model-based controllers. ISA Transactions, 47:395–400, 2008.[21] M. Yousefi, L.D. Rippon, and M.G. Forbes. Moving-horizon predictive input designfor closed-loop identification. 9th International Symposium on Advanced Control ofChemical Processes, pages 135–140, 2015.84
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A novel algorithm for model plant mismatch detection...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A novel algorithm for model plant mismatch detection for model predictive controllers Tsai, Yiting 2016
pdf
Page Metadata
Item Metadata
Title | A novel algorithm for model plant mismatch detection for model predictive controllers |
Creator |
Tsai, Yiting |
Publisher | University of British Columbia |
Date Issued | 2016 |
Description | For Model Predictive Controlled (MPC) plants, the quality of the plant model determines the quality of performance of the controller. Model Plant Mismatch (MPM), the discrepancies between the plant model and actual plant transfer matrix, can both improve or degrade controller “performance”, depending on the context in which “performance” is defined. Instead of using simple, “yes-no” type of performance metrics to diagnose whether MPM is present, this thesis achieves the further goal of pinpointing the exact plant inputs causing MPM. The detection method consists of a two-step identification procedure. The first experiment identifies the MPM-affected rows in the plant matrix. The second experiment pinpoints the exact inputs causing said MPM, and these MPM-affected elements are further ascertained using a hypothesis test. Finally, using input design, an estimated optimal excitation sequence is generated based on the available closed-loop data, which helps the engineer determine whether the original excitation was sufficient for the aforementioned sys-ID experiments. The most important underlying assumptions are the linear time-invariance of the plant and noise dynamics. The proposed algorithm is exercised on artificial 3x3 and 5x5 plants, then on real data from an industrial lime slaker suffering from sparse MPM, to demonstrate its ability of accurately pinpointing MPM-affected elements. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2016-04-26 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0300243 |
URI | http://hdl.handle.net/2429/57859 |
Degree |
Master of Applied Science - MASc |
Program |
Chemical and Biological Engineering |
Affiliation |
Applied Science, Faculty of Chemical and Biological Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2016-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2016_may_tsai_yiting.pdf [ 2.33MB ]
- Metadata
- JSON: 24-1.0300243.json
- JSON-LD: 24-1.0300243-ld.json
- RDF/XML (Pretty): 24-1.0300243-rdf.xml
- RDF/JSON: 24-1.0300243-rdf.json
- Turtle: 24-1.0300243-turtle.txt
- N-Triples: 24-1.0300243-rdf-ntriples.txt
- Original Record: 24-1.0300243-source.json
- Full Text
- 24-1.0300243-fulltext.txt
- Citation
- 24-1.0300243.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0300243/manifest