Stabilizing

We demonstate stabilization of a computational fluid dynamics model of an unstable system. The unstable heating of a two-dimensional plate is used as a case study. Active control is introduced by cooling parts of the boundaries of the plate. The high order of the original model is reduced by proper orthogonal decomposition, giving an unstable reduced order model with a state space structure convenient for controller design. A stabilizing controller based on pole placement is designed for the reduced order model and integral action is included to enhance performance. The controller is then applied to the full model, where it is shown through simulations to stabilize the system. The demonstrated procedure makes it possible to analyze stability properties and design control systems for a class of systems that would otherwise be very computationally demanding.


Introduction
Control of ‡uid ‡ow is a research area that has gained an increasing amount of interest over the last decade, with a recent contribution in [1].Both internal ‡ow ( ‡ow inside bounded regions, such as within engines or turbomachinery) and external ‡ow ( ‡ow around vehicles or bodies) have recently been subject to an increasing amount of interest from the control community.A nice review of advances made in the intersection between control theory and ‡uid mechanics can be found in [2].Distributed systems such as ‡ow dynamics are modeled mathematically by a set of partial differential equations (PDEs).In order to simulate such systems on a computer, the PDEs are usually discretized using a set of tools known as computational ‡uid dynamics (CFD) ( [3], [4], [5]).Although often acceptable for simulation purposes, the resulting discretized models are often of such high order that it may be infeasible to use them in controller design.For instance, numerical simulation of the full three-dimensional Navier-Stokes equations is still too costly for the purpose of optimization and control of unsteady ‡ows [6].
There is, however, a potential for using model order reduction tools to reduce the complexity of the models.Model reduction is a powerful tool, and a recent survey of di¤erent techniques may be found in [7].Although many of them have been successfully used within the control community, methods such as Hankel model reduction and balanced truncation are too computationally demanding for systems of such high order as those encountered in many CFD applications.
Proper orthogonal decomposition (POD) is a popular technique that is still applicable even for very high-order systems, using the method of snapshots that will be further presented in section 2. Nice examples of this approach have been shown e.g. in [6] and [8], in active control of ‡uid ‡ow over a backward-facing step and forward-facing step, respectively, using …nite element discretization of the Navier-Stokes equations.In [9] and [10] it was shown that the reduced order models can have a state space structure, which is convenient for controller design.In this paper we extend the analysis of the unsteady two-dimensional heat equation studied in [10] to also include a heat source that renders the plate model unstable.
The paper is organized as follows.In section 2 the model order reduction technique is presented.The full CFD model of the unstable system is derived in section 3.1, and the reduced order model is found in section 3.2.A controller is designed in section 4 and validated in section 5, before concluding remarks are o¤ered in section 6.
First introduced independently by Karhunen [11] and Loève [12], POD is sometimes called the Karhunen-Loève expansion.When …rst introduced in the context of ‡uid mechanics in [13], it was used to study turbulent ‡ows.Applicable even for complex non-linear problems, POD makes it possible to keep the essential dynamics of the original model in order to perform the controller design.We will now outline the procedure.Consider the general transient model of a physical system @x @t = D (x) ; where x 2 R n is the state, is a spatial coordinate and D (x) is some nonlinear partial di¤erential operator in the spatial variables [14].By assuming that the state x can be represented as a linear combination of r basis functions, where x r 2 R r is an approximation of the state, a reduced order model of ( 1) can be derived.The matrix r 2 R n r is an orthogonal projection matrix (i.e.T r r = I) whose columns are orthonormal basis vectors i that span the 'space'of spatial dynamics of the state, and the i 's are time varying expansion coe¢ cients (also called modal amplitudes).The POD o¤ers a way to …nd a reduced set of basis functions i , based on an ensemble of data of the state trajectories x (t), using the method of snapshots [15].The snapshots may be taken from physical experiments or from computer (CFD) simulations.Since the procedure is purely data-driven, i.e. the basis functions are found without considering the governing equations, …nding the snapshot data is critical, since the reduced-order model will capture only the dynamics present in the data.A suitable input should therefore be used to excite the system, so that the desired characteristics are present in the data.How to best choose suitable inputs is discussed in [10].The snapshots are collected in a matrix X , X = x 1 ; x 2 : : : ; x N = [x (t 1 ) ; x (t 2 ) ; : : : ; x (t N )] ; (3) where the columns fX ; j g N j=1 can be thought of as the spatial coordinate vectors of the system at time step t j .The rows fX i; g n i=1 describe the time trajectories of the system evaluated at di¤erent locations in the spatial domain [16].The POD basis vectors can be found by choosing the orthonormal set that maximizes the cost [17] where ( ; ) is the scalar product operator and h i is a time averaging operator.Equivalently, the basis functions can be found by performing singular value decomposition (SVD) of the snapshot matrix X , which is the procedure followed in this paper.A rigorous treatment of the connection between SVD and POD can be found in [18].By SVD, the snapshot matrix is decomposed as where the columns of = [ 1 ; : : : n ] form the optimal orthogonal basis for the space spanned by X , and are unitary matrices (i.e. 1 = T ; 1 = T ) and P is a diagonal matrix with the singular values of X on the diagonal.The r most signi…cant basis functions are associated with the r largest singular values of X .If the singular values fall of rapidly, a reduced order model may be constructed that captures the most salient characteristics of the original system.The time varying expansion coe¢ cients i in equation (2) are determined through Galerkin projection, a well-known method for generating a system of ordinary di¤erential equations from a system of partial di¤erential equations [19], [20].To ensure that the original model is approximated as accurately as possible, i are found by requiring that the projection of the truncated solutions onto the basis is zero, i ; @x r @t D (x r ) = 0; i = f1; : : : rg : By posing this constraint on i (t) for a given set of basis functions i , the expansion coe¢ cients i can be found by inserting (2) into (5) to yield r initial value problems for : If the initial values of equation ( 1) are given as x (0; ) = f ( ), where is a spatial variable, the initial value for is found from Thus the problem of solving a possibly complex PDE is reduced to that of solving a set of ordinary di¤erential equations.Eventually, pursuing a reduced order system on the regular linear state space form, the model reduction objective is to …nd a model where where r n, that approximates the system (1).
3 Case study: Heated plate

Full CFD model
To demonstrate how an unstable system can be stabilized using POD and feedback control, we study heat conduction in a plate.The plate is 1 m 1 m, 1 cm thick, de…ning the twodimensional 1  1 A two-dimensional analysis is valid since the plate is so thin.
the plate is governed by the unsteady linear twodimensional heat equation where and c p are the density and speci…c heat capacity of the plate, respectively, and k is the thermal conductivity, that is assumed to be uniform over the computational domain and independent of the temperature.Note that x now and in the following denotes a spatial coordinate and no longer the state variable.The term S , S c + S T describes heat sinks and sources.In the present problem, convective heat transfer to the surroundings gives rise to a term where h is the convective heat transfer coe¢ cient, A is the heat transfer area of the surface and T 1 is the ambient temperature.Expression ( 10) is a source term whenever T > T 1 , and a sink term whenever T < T 1 .Due to electric current, the plate is subject to an internal temperaturedependent heat source where k 1 > 0, at all points except from the boundary.Intuitively, this positive feedback from the temperature to the source may lead to a physically unstable system if the convective heat loss to the surroundings is not large enough.An increase in temperature will then lead to a stronger source, which again increases the temperature, and so on.Discretizing the governing equation by the …nite volume method, ( 9) is integrated over each control volume CV and over the time interval from t to t + t, to obtain [5] where the order of integration has been changed for the rate of change term.Using the unconditionally (numerically) stable backward Euler (fully implicit) temporal discretization and M grid points over the spatial domain , the system (9) can be written as a system of M equations on the form ) where the a's are coe¢ cients and T P is the temperature at the grid point (point P ) under consideration at time step k + 1. S u and S P arise from discretizing the source term S as V S = S u +S P T P .Using the convenient compass notation, T W , T E , T S and T N are the temperatures at the west, east, south and north neighbor grid points, respectively, at time step k + 1. T 0 P is the temperature at grid point P at time step k.Collecting the temperature at all grid points in a row vector T (k) 2 R M leads to a discrete linear system on descriptor form where E 2 R M M is a penta-diagonal matrix containing the coe¢ cients a p , a W , a E , a S and a N and A 2 R M M is a diagonal matrix with a 0 P on the main diagonal.B 2 R M m contains the contributions from the inputs, while the constant source terms give rise to a constant term V 2 R M .Linear systems such as (12) can be solved either by direct methods, such as matrix inversion or Gauss elimination, or by iterative methods, such as successive over-relaxation or more sophisticated solvers.In the case of non-linear systems, for instance when the thermal conductivity k is dependent on the temperature or the source term is nonlinear, iteration is required.In the CFD community, iterative methods are often preferred rather than direct methods even when the equation system is linear, mainly due to signi…cantly smaller storage requirements.In a domain comprising M grid points, a direct method will require the simultaneous storage of M 2 coe¢ cients, while an iterative solver will only require the storage of the non-zero coe¢ cients.Given the sparsity of the co-e¢ cient matrices, the pro…t is considerable.To validate that the plate model is unstable, we rearrange system (12) into regular linear state space form by inverting E and multiplying throughout (12) to obtain where The matrix E is indeed invertible in this case, since it is strictly diagonally dominant, and consequently non-singular.This is not necessarily possible in general, but it allows us to establish that the system matrix A has an unstable eigenvalue outside the unit circle, j (A) j max = 1:00001 > 1: (14)

Reduced order model
The PDE ( 9) is discretized using 50 grid points in both the x-and y-direction.This gives in total 2500 states in the full model (13).To construct a model of reduced order, ( 13) is simulated for 600 time steps, thus forming the matrix of snapshots X .During this simulation the inputs are varied to excite as much of the dynamics of the system as possible.SVD of the snapshot matrix is performed, and the …rst four columns of , corresponding to the four largest singular values are extracted to form the POD basis r , as depicted in …gure 2. As can be seen from the …gure the singular values fall o¤ quite rapidly, and many of the singular values are zero, indicating that the basis functions corresponding to those singular values can be omitted without loss of information.With  the …rst four basis functions, the heuristic crite-rion indicating a fairly accurate reduced order model [9].The criterion gives an indication of the part of the energy that is conserved in the reduced order model.Moreover, if the reduced order state has dimension four, then the number of states in the reduced order model is equal to the number of inputs.Consequently, the reduced order model is fully actuated, which might be favorable when tracking a reference pro…le for the complete state.The basis r de…nes the linear dimensionalityreducing mapping H : R M 7 !R r , such that = H (T ) = T r T , and the expansion G : R r 7 !R M such that T = G ( ) = r .To obtain the reduced order model on the form (8), the procedure presented in [9] is followed.By applying the inner product criterion ( 6) and inserting T = into (12), the reduced order model can be obtained: ) De…ning E r , T E r allows us to write where E r is invertible since E, T and r are all nonsingular.This yields the reduced-order model on discrete state space form where and Ĉ 2 R p r .In this example, r = m = p = 4, and Ĉ is chosen to be I 4 4 .In an actual implementation, i should be reconstructed from measurements on the real state T (t; x; y) through an observer, rather than measuring the whole reduced state, as done in [9].Here we merely assume that all states i are available for feedback.The reduced order model is unstable, since that is, outside the unit circle.In fact, the unstable eigenvalue of the reduced-order model is larger than the unstable eigenvalue of the full model.Note that the general POD procedure does not automatically preserve stability properties during the model reduction process.Criteria for preserving stability properties are derived in [21].

Controller design
Feedback control is performed by use of heat ‡ux actuators on parts of the boundary of the domain.
The actuators are placed at the bold lines in …gure 1.The control objective is to reach a constant temperature reference T d while at the same time rejecting disturbances.The reference temperature is set to be a uniform temperature of 300 K. Since the full model is too large for controller design the reduced order model is analyzed instead.The reduced order state reference d is found as d = T T d .Given the unstable reduced order model (18), the control objective is to stabilize the system around the reference temperature.
De…ning the tracking error e (k) , d y (k), and recalling that y (k) = Ĉ (k) = (k), the control input is chosen as where K is chosen such that the eigenvalues of the closed-loop system matrix Â BK have magnitude less than one.Taking into consideration that the reduced-order model is merely an approximation, the controller should include integral action in order to minimize the steady-state error.To do this in a straightforward way, we de…ne the augmented state giving an augmented state space model where and u (k) = u (k) u (k 1).In this augmented state space, integral action is built-in, and the input increment u (k) is found as where K is the feedback gain matrix found above.
Without control the temperature of the plate is strictly increasing.After 20 minutes the highest plate temperature is 368 K.If the simulation is run for a longer period of time the temperature continues to increase, illustrating the instability of the system.Now, the full CFD-model is simulated with the controller designed for the reduced order model in section 4. Initially, the plate temperature is equal to the ambient temperature at 293 K.At t = 0 the inner source and controller are switched on, and the system is simulated until steady-state occurs, after approximately 20 minutes.The largest steady-state error is approximately 4 K, as shown in …gure 3.  It is seen that although the original CFD model is symmetric, the controller based on the reducedorder model does not manage to exploit this symmetry, since the symmetry is not preserved in the model reduction scheme.The simulation of the full closed-loop CFD model illustrates that the plant has been stabilized.The temperature distribution at steady-state is shown in …gure 4.

Concluding Remarks
In this paper we have demonstrated, using a case study, that a CFD-model of an unstable system can be stabilized through model reduction and a controller designed for the reduced order model.This makes it possible to study stability proper-  ties of systems that would otherwise be very computationally demanding.Further and ongoing work include stability analysis of reduced-order models of more general PDEmodels, and development of model reduction methodology.The 2D-plate model of the current work is viewed as a step on the way in this work.

Figure 1 :
Figure 1: Sketch of plate with actuators on boundaries (bold lines).

Figure 2 :
Figure 2: Singular values of the snapshot matrix.The 's indicate singular values corresponding to the extracted basis functions.Note that the ordinate axis is logarithmic.From singular value 10-600, every tenth singular value is plotted, for clarity.