Thermal vacuum chamber control optimization

TREBALL DE FI DE CARRERA TÍTOL DEL TFC: Thermal Vacuum Chamber Control Optimization TITULACIÓ: Enginyeria Tècnica Aeronàutica, especialitat Aeronavegació AUTOR: Dídac Jabaloy Colominas DIRECTOR: Pr. Franco Bernelli Zazzera DATA: 27/02/2012 Títol: Thermal Vacuum Chamber Control Optimization Autor: Dídac Jabaloy Colominas Director: Pr. Franco Bernelli Zazzera Data: 27/02/2012 Resum En un laboratori del departament d’Enginyeria Aeroespacial de la Universitat Politècnica de Milano (Itàlia) hi ha una càmera de buit tèrmic. L’objectiu d’aquest tipus de màquines és testejar compo nents espacials. La màquina està sent millorada per projectes d’estudiants. Aquesta tesis té l’objectiu d’incloure òptimes capacitats de control per l’operació de la màquina. La forma típica de treballar de la càmera de buit tèrmic és la següent. Es posa un objecte dins seu, aleshores la càmera es segella i es treu tot l’aire. Des d’aquest punt, l’objecte pot se refredat o escalfat. El procés d’escalfament està controlat per un contr olador PID mentre el de refredament ho està per un sistema de control estil onoff. Per aquest motiu, l’estudi de l’escalfament i el refredament de la càmera és completament divers. Com l’escalfament té més flexibil itat per a ser optimitzat, aquest projecte te l’objectiu d’estudiar un possible canvi del controlador PID a un nou microcontrolador que millori les característiques del sistema en els termes de consum d’energia. L’idea es dissenyar un mètode de treball més suau i estable que redueixi el consum elèctric. A més, amb la potència de càlcul dels microcontroladors actuals, és possible dissenyar un algoritme que canviï els seus paràmetres d’acord amb les condicions de treball de la màquina. És a dir, que quan la Càmera comença el procés, s’analitza el seu funcionament i s’escullen els millors paràmetres per treballar en aquell moment. S’ha de tenir en compte que les característiques de l’escalfament de la Càmera canviaran cada vegada a causa de les propietats tèrmiques de l’objecte que és dins seu. Per dur a terme l’estudi s’ha utilitzat la teoria de control de màquines físiques mitjançant la predicció de la sortida amb sistemes d’equacions diferencials. Aquesta teoria matemàtica està basada en el fet que, per q ualsevol sistema físic, és possible predir la seva resposta (sortida o output) a qualsevol estimulació (entrada o input). Tot el software i els gràfics utilitzats en aquesta tesis han estat escrits utilitzant l’entorn de càlcul numèric MATLAB. El futur co di a compilar en el microcontrolador serà una compilació dels codis utilitzats en aquest treball. Title: Thermal Vacuum Chamber Control Optimization Author: Dídac Jabaloy Colominas Director: Pr. Franco Bernelli Zazzera Date: 27/02/2012 Overview In a laboratory of the Aerospace Engineering department of the Polytechnic University of Milano (Italy) there is a thermal vacuum chamber. The objective of this type of machines is to test space components. The machine is being improved by student’s projects. This thesis has the objective of including optimal control capabilities for machine’s operation. The typical way of working of the thermal vacuum chamber is the following. An object is put inside the chamber. Then, the cham ber is sealed and all air is realised. From this point, the object could be cooled or heated. The heat process is controlled by a PID controller while the cooling is controlled by an on-off control. For this reason, the study of chamber’s heating and cool ing is completely different. How the heat has more flexibility to be optimized, this project has the aim of studying a possible change of the PID controller for a new microcontroller that would improve the performances of the system in terms of energy consumption. The idea is to design a smoother and more stable working method that reduces the electricity consumption. Furthermore, with the nowadays microcontrollers calculus power, it is possible to design an algorithm that changes its parameters according to machine’s working conditions. That is to say that, when the Chamber starts the process, are analysed its performances and are chosen the best parameters to work in that moment. Remember that the Chamber heating performances will change every time because of the thermal characteristics of the object that is inside. To carry out the study it has been used the theory of control of physical machines by output prediction with difference equations systems. This mathematical theory is based on the fact that, f or any physical system, it is possible to predict its response (output) to any stimulation (input). All software and graphics used in this thesis has been written using the numerical computing environment MATLAB. The future code to be compiled in the microcontroller will be a compilation of the codes used in this work. CONTENTS Introduction ....................................................................................................... 11 1. Thermal Vacuum Chamber ....................................................................... 13 1.1 Vacuum Circuit .................................................................................... 13 1.2 Heating Circuit ..................................................................................... 14 1.3 Cooling Circuit ..................................................................................... 14 1.4 Sensors and Data Acquisition ............................................................. 15 1.5 Control System .................................................................................... 15 2. System Identification ................................................................................. 17 2.1 Finite-difference model equation ......................................................... 17 2.2 Getting alpha & beta ............................................................................ 18 2.3 Alpha and Beta of the Thermal Vacuum Chamber .............................. 19 3. Multistep Output Prediction ....................................................................... 21 3.1 Multistep Single Function .................................................................... 21 3.2 Multistep Matrix Function .................................................................... 22 3.3 Ideal system order and prediction horizon ........................................... 23 4. Multistep Output Prediction Application ..................................................... 25 4.1 Calculating the prediction error ........................................................... 25 4.2 Prediction Example with single input and output ................................. 26 4.3 Prediction Example with Multiple Inputs and Outputs.......................... 29 4.4 Prediction Application on the Thermal Vacuum Chamber ................... 31 5. Predictive Control Algorithm ...................................................................... 35 5.1 Generalized Predictive Control Theory ................................................ 35 5.2 GPC Parameters Theoretical Effect .................................................... 36 5.3 Energy Consumption Calculus ............................................................ 37 6. Control algorithm application ..................................................................... 39 6.1 Control Example with single Input and Output .................................... 39 6.2 Control Example with Multiple Inputs and Outputs .............................. 41 6.3 Control Application on the Thermal Vacuum Chamber ....................... 44 8 Thermal Vacuum Chamber Control Optimization 7. Conclusions ............................................................................................... 47 Bibliography ...................................................................................................... 49 ANNEXES ........................................................................................................ 51 A. Pictures about the Chamber ................................................................ 53 B. Alpha and Beta Calculation ................................................................. 55 C. Multistep Matrix Generator with Single input/output ............................ 57 D. Multistep Matrix Generator with Multiple input/output .......................... 59 E. Example with Multiple Inputs and Outputs .......................................... 63 F. Generalised Predictive Control Algorithm ............................................ 65 LIST OF FIGURES Figure 0.1 PID Controller Output Voltage ......................................................... 11 Figure 1.1 The Thermal Vacuum Chamber ...................................................... 13 Figure 1.2 Turbomolecular Pump ..................................................................... 14 Figure 1.3 Inside the Chamber ......................................................................... 15 Figure 1.4 Control Panel .................................................................................. 16 Figure 2.1 Output Example with the Chamber Empty ...................................... 20 Figure 4.1 Prediction operation example .......................................................... 25 Figure 4.2 Scheme of the first example ............................................................ 26 Figure 4.3 Output prediction difference with p=1 and p=2. .............................. 27 Figure 4.4 Error in function of the system order. .............................................. 28 Figure 4.5 Error in function of prediction horizon. ............................................. 28 Figure 4.6 Scheme of the example 2. ............................................................... 29 Figure 4.7 Real output vs. predicted with p=2 and s=5. ................................... 30 Figure 4.8 Error in function of the system order. .............................................. 30 Figure 4.9 Error in function of prediction horizon. ............................................. 31 Figure 4.10 Surface 3D of P, S and Error ......................................................... 32 Figure 4.11 Combination of P-S that satisfies the maximum error. .................. 32 Figure 4.12 Real output vs. predicted in a test from the Chamber. .................. 33 Figure 5.1 Energy Calculation .......................................................................... 38 Figure 6.1 Example 1 controlled with p= hc=2 and w=10. .............................. 39 Figure 6.2 Error in function of weight ................................................................ 40 Figure 6.3 Energy in function of weight ............................................................ 41 Figure 6.4 Input with p=2, hc=2, w=10-2 ........................................................... 42 Figure 6.5 Output with p=2, hc=2, w=10-2 ......................................................... 42 Figure 6.6 Input with p=2, hc=2, w=0.5 ............................................................. 43 Figure 6.7 Output with p=2, hc=2, w=0.5 .......................................................... 43 Figure 6.8 Old versus controlled Chamber’s input ........................................... 44 Figure 6.9 Old versus controlled Chamber’s output ......................................... 45 LIST OF TABLES Table 2.1 Finite-difference model equation’s terms. ......................................... 17 Table 2.2 Chamber’s alphas and betas with an order model of 5. .................... 20 Table 5.1 Change in the control depending on its parameters. ........................ 37 16 Thermal Vacuum Chamber Control Optimization The on-off control id used for the cooling circuit because it is difficult to know the exact Freon power absorbed in the serpentine. On the other hand, for the heat circuit, the PID control is used, and the controllers receive a voltage that is proportional to the difference between the reference temperature and the real temperature inside the chamber and its derivative. Figure 1.4 Control Panel To see more pictures about the Thermal Vacuum Chamber go to the Annex A. Pictures about the Chamber. System Identification 17 2. SYSTEM IDENTIFICATION In this context, identification is the process of developing a mathematical model made from experimental data to apply in a physical system. Are considered a physical system all structures or machines from which a logical output is obtained after applying an input. In this essay, the system is the thermal vacuum chamber. The input is the voltage imposed by the heat controller, and goes to the lamps and resistances controllers, and the output is the temperature inside the chamber. First, the relationship between the input and output data will be obtained in the form of a time-domain finite difference model. This model is known as the autoregressive model with exogenous input (ARX). Next, in this chapter, the coefficients of the finite-difference model from input and output data will be computed. 2.1 Finite-difference model equation The relation between the input and the output in a physical system is usually described by a finite-difference model equation. This equation and his coefficients are commonly (to find more information about them, see the bibliography [1]): y(k)=−α1· y(k−1)−α2· y(k−2)−⋯−αp· y(k−p)+ +β0· u(k)+β1· u(k−1)+⋯+βp· u(k−p). (2.1) Table 2.1 Finite-difference model equation’s terms. Coefficient Description u(k) Input time history (where r = nº of inputs) y(k) Output time history (where m = nº of outputs) p Integer indicating the order of the model α1… αp Matrixes [m×m] associated with the [m×l] output vector 𝐸 at time steps (k−1),…,(k−p) respectively 18 Thermal Vacuum Chamber Control Optimization β0,β1… βp Matrixes [m×r] associated with the [r×l] input vector 𝐸 at time steps (k), (k−1),…,(k−p) respectively The equation uses a 𝑝 number of past outputs and input measurements y(k−1), … , y(k−p)and u(k−1), … , u(k−p), and the current input measurement u(k) to get the output in progress y(k). Note that equation (2.1) refers to the arbitrary 𝑘 time step. This process is a forward-in-time prediction and is easy to do, knowing the alpha and beta coefficients. Unfortunately, these coefficients are usually unknown, so they must be estimated. Another important topic of this equation is about the single or multiple input/output. All the theory exposed in this thesis works from 1 to any number of inputs/outputs because the target is that the code could work correctly for any system. Nonetheless, the vacuum chamber’s equation is modeled with one input and one output. It is important to understand well these concepts and their size since they are the baseline of all following mathematical development. 2.2 Getting alpha & beta With the background of inputs and outputs history, it is possible to obtain the alpha and beta coefficients. This mathematical process is complex and extremely hard without computer calculus. Firstly is defined the matrix v(k) and its length between ⟦ ⟧: v(k)=�y(k) u(k)�,(k = 1, 2, … , l); ⟦(m + r)× l⟧ (2.2) where, l is the length of data. If all outputs in equation (2.1) are moved to the left side of the equal, with k from 1 to l, the following matrix equation is obtained: [y0 y]=θ[V0 V] (2.3) The length of [V0 V] is ⟦[p · (m + r) + r] × l⟧ . The (2.3) factors are: y0= [y(1) y(2)… y(p)] (2.4) y = [y(p + 1) y(p + 2)… y(l)] (2.5) System Identification 19 θ= [β0 (−α1 β1)… �−αp−1 βp−1� �−αp βp�] (2.6) V0= ⎣ ⎢ ⎢ ⎢ ⎡ u(1) u(2) ⋯u(p) v(0) v(1) ⋯v(p −1) ⋮ ⋮ ⋱ ⋮ v(2 −p) v(3 −p) ⋯v(1) v(1 −p) v(2 −p) ⋯v(0) ⎦ ⎥ ⎥ ⎥ ⎤ (2.7) V = ⎣ ⎢ ⎢ ⎢ ⎡ u(p + 1) u(p + 2) ⋯u(l) v(p) v(p + 1) ⋯v(l −1) ⋮ ⋮ ⋱ ⋮ v(2) v(3) ⋯v(l −p + 1) v(1) v(2) ⋯v(l −p) ⎦ ⎥ ⎥ ⎥ ⎤ (2.8) Obviously, we are interested in finding θ. Thus, the pseudoinverse (also called generalized inverse or Moore–Penrose pseudoinverse) is computed in order to get: θ=[y0 y]·[V0 V]†P (2.9) A good way for simplify this equation is to work as if the system in study was with initial transient response. That is to say, until the time step 1 there are not input neither output different from 0. By this way it is possible to eliminate V0 and y0. The better point is that this new equation will be applicable in the case without transient response because it founds some coefficients that describes a system and does not matter if a short part (time steps since 0 till p) is not included. So the final equation to get alpha and beta is: θ= y · V†P (2.10) 2.3 Alpha and Beta of the Thermal Vacuum Chamber It is created a code that founds the alpha and beta in the chamber using a value of p chosen by the user (See annex B. Alpha and Beta Calculation). The rationales for choosing this value will be discussed later. With the machine empty, and using all components like in a real test, was made a test in order to get experimental data. This test consisted in an input impulse, 20 Thermal Vacuum Chamber Control Optimization of maximum voltage, of 5 minutes duration with some more minutes without any input to leave the Chamber be cooler. The output (the temperature of the Chamber) was measured only in one sensor, whose results are the following: Figure 2.1 Output Example with the Chamber Empty With the impulse input and this output is applied the code to get the p terms with an order model of 5. The coefficients obtained are in the next table: Table 2.2 Chamber’s alphas and betas with an order model of 5. Subindex α β 0 - 0.0084 1 -0.4741 -0.0050 2 -0.3121 -0.0003 3 -0.1843 -0.0039 4 -0.0825 0.0010 5 0.0529 0.0108 Multistep Output Prediction 21 3. MULTISTEP OUTPUT PREDICTION In the previous chapter, it is explained how to get the present output by the present input and past inputs and outputs (formula (2.1)). If the wanted output is not the next step in the time domain, that formula needs some improvements. This chapter illustrates how to find the future output of a physical system at any specific time step. 3.1 Multistep Single Function The basic idea is the same, rests simply on the relationship between the current output and the previous inputs and outputs signals. For example, the output on the step k+1 is only the equation (2.1) but shifting the index one step ahead: y(k + 1)=−α1· y(k)−α2· y(k−1)−⋯−αp· y(k−p + 1)+ +β0· u(k + 1)+β1· u(k)+⋯+βp· u(k−p + 1). (3.1) Note that in this equation is required the output signal at one time step aback, so actually the prediction is only for one step yet. Substituting y(k) from equation (2.1) and also shifting one step the alpha and beta coefficients (the x(1) says that are for one step ahead) with a method showed later: y(k + 1)=−α1(1)· y(k−1)−α2(1)· y(k−2)−⋯−αp(1)· y(k−p)+ +β0· u(k + 1)+β0(1)· u(k)+β1(1)· u(k−1)+⋯+βp(1)· u(k−p). (3.2) Now is developed an equation where the output is obtained on the step k+1 without knowing it at k . Notice that the future input signals must be given a priori. If the process to arrive at (3.2) is applied many times, the output measurement at time step k+j will be expressed as: y(k + j)=−α1(j)· y(k−1)−α2(j)· y(k−2)−⋯−αp(j)· y(k−p)+ +β0· u(k + j)+⋯+β0(1)· u(k + j −1)+⋯+β0(j)· u(k)+ +β1(j)· u(k−1)+β2(j)· u(k−2)+⋯+βp(j)· u(k−p). (3.3) The coefficients must be shifted for an amount of steps equal to the prediction step to be predicted. To shift them, the next functions have to be performed: 22 Thermal Vacuum Chamber Control Optimization β0(0)=β0 β0(j)=βj−� αiβ0(j−i) ji=1 for j = 1, 2, … , p. β0(j)=−� αiβ0(j−i) ji=1 for j = p + 1, … , ∞. (3.4) αi(0)=αi αi(j)=αi+1 (j−1)−α1(j−1)·αi for i = 1, 2, … , p −1. αp(j)=−α1(j−1)·αp (3.5) βi(0)=βi βi(j)=βi+1 (j−1)−α1(j−1)·βi for i = 1, 2, … , p −1. βp(j)=−α1(j−1)·βp (3.6) 3.2 Multistep Matrix Function The solution of this equation (3.3) is the output at a specific point in the future. More interesting would be to find many points at same time in order to obtain the output’s evolution, a prediction vector or matrix of outputs. In this way, if the index j is let to be j=1, 2, …, s-1, the equation (3.3) will produce the matrix equation: ys(k)= T · us(k)+ B · up(k−p)−A · yp(k−p) (3.7) Whose constant matrices T, B, A and the other components are computed by (their length is between ⟦ ⟧): T = ⎣ ⎢ ⎢ ⎢ ⎡ β00 0 0 β0(1)β00 0 ⋮ ⋮ ⋱ 0 β0(s−1)β0(s−2)⋯ β0 ⎦ ⎥ ⎥ ⎥ ⎤ ; ⟦𝑠·𝑚×𝑠·𝐸⟧ (3.8) Multistep Output Prediction 23 B = ⎣ ⎢ ⎢ ⎢ ⎡ βpβp−1 ⋯ β1 βp(1)βp−1 (1)⋯ β1(1) ⋮ ⋮ ⋱ ⋮ βp(s−1)βp−1 (s−1)⋯ β1(s−1) ⎦ ⎥ ⎥ ⎥ ⎤ ; ⟦s · m × p · r⟧ (3.9) A = ⎣ ⎢ ⎢ ⎢ ⎡ αpαp−1 ⋯ α1 αp(1)αp−1 (1)⋯ α1(1) ⋮ ⋮ ⋱ ⋮ αp(s−1)αp−1 (s−1)⋯ α1(s−1) ⎦ ⎥ ⎥ ⎥ ⎤ ; ⟦𝑠·𝑚×𝑝·𝑚⟧ (3.10) us(k) = �u(k) u(k + 1) ⋮ u(k + s −1)� (3.11) up(k −p) = �u(k −p) u(k −p + 1) ⋮ u(k −1) � (3.12) The expressions for the outputs, ys(k) and yp(k−p), are the same as for the inputs (equations 3.11 and 3.12) but changing u for y. The s is an integer related to the prediction horizon: hs=s-1. The prediction horizon is the last point in the future that is being predicted. The study of this parameter is developed in the following section. 3.3 Ideal system order and prediction horizon In a differential equation, the order is the maximum derivative grade of the equation. Logically, the higher is the order the higher is the complexity that can be the system that describes the equation. With real physical systems, there are so many variables that influences in the responses that is not possible to say exactly the differential equation order. Furthermore it is important not to forget that some quantity of noise always affects real situations. This noise will deform randomly the real output, adding more complexity to understand the relation between input and output. For these reasons, it is not unreasonable to think that the best option would be to use the maximum system order possible. This premise is correct a priori, but unacceptable in a real project. Choosing a very high system order, the 24 Thermal Vacuum Chamber Control Optimization computational time for output prediction would not compensate the error reduction. On the other hand, the parameter s (prediction horizon) works in the contrary way, it has an accumulative error. The longest period you want to predict, the biggest error you get in comparison with the real output. Another study has been performed to find the optimum prediction horizon (hs=s-1). Multistep Output Prediction Application 25 4. MULTISTEP OUTPUT PREDICTION APPLICATION How it is seen in the previous chapter, for the multistep output prediction there are two variables to be chosen: the system order (p) and the prediction horizon (hs). In order to choose the better combination of them, it is essential to study experimentally the changes that they produce. In this chapter, first will be presented different examples to test that the implemented software works correctly; that is, if parameter p is increased the prediction error respect the real output will decrease and just the contrary with the parameter s. Secondly will be applied all the contents of chapters 3 and 4 in the thermal vacuum chamber to get the optimum p and s for this thesis. 4.1 Calculating the prediction error To understand better the way to get the prediction error, an explanation of how the program works with the prediction horizon follows: - Firstly, the program wait that the system produces p number of outputs because this is the minimum number to calculate the parameters alpha and beta that describes the system. - Next is made the first prediction. From the point p, are calculated s-1 number of outputs. The maximum error is in the last step calculated, so the interest is there. From the second prediction until the end of data, only is taken the last step prediction every time. The next figure illustrates all the process. Figure 4.1 Prediction operation example 32 Thermal Vacuum Chamber Control Optimization Figure 4.10 Surface 3D of P, S and Error. Next figure shows all combinations of p and s whose prediction error is under 0.1ºC, so always should be taken some of these points to work with the Chamber: Figure 4.11 Combination of P-S that satisfies the maximum error. Multistep Output Prediction Application 33 If for example, with some real data from the Chamber in a test which is reach a temperature of 40ºC is applied an output prediction with p=11 and s=6. The result shows that the predicted output is almost exactly like the real output: Figure 4.12 Real output vs. predicted in a test from the Chamber. Predictive Control Algorithm 35 5. PREDICTIVE CONTROL ALGORITHM After learning the formula that give us the outputs of a system at any series of points in the future, it is the moment to think about the input. Until now, the input, also called control, is given a priori to estimate the output with it. In this chapter again is predicted the output, but also the needed input to get a desired output is calculated. A generalized predictive control algorithm is used to carry out this process. To find more information about this algorithm and its formulas, see the bibliography [2]. 5.1 Generalized Predictive Control Theory The Generalized Predictive Control (GPC) algorithm is based on the multistep output prediction. It is computed with the matrixes formed from the stepresponse time history in conjunction with a cost function. The control input is obtained by minimization of the cost function. For these reasons there are some rules about the future control inputs: - The control is applied only till the control horizon (hc), which must be equal to or less than the prediction horizon (hp). After the control horizon, the control is assumed to be zero. Actually, in this thesis, is assumed that: hc=hp. - In addition, the function has a control penalty, a “weight”, that give us the possibility of change some control properties and stabilize the closedloop system. This is because it gives more importance to error minimization or control minimization depending on its value.. The cost function to be minimized in the GPC algorithm is: J(k)=12 �·{[ys(k)−y�s(k)]T·[ys(k)−ys(k)]+ usT(k)·λ· ys(k)} (5.1) Where ys(k) is the predicted output vector seen in last chapter, while y�s(k) is the desired output vector. The weight or control penalty is λ and it is a positive scalar. Minimizing the equation with respect to us(k) yields: ∂J(k) ∂us(𝑘)=�∂ys(k) ∂us�T·[ys(k)−ys �(k)] +λ· us(k)≡0 (5.2) Note from equation (3.7) that: 36 Thermal Vacuum Chamber Control Optimization ∂ys(k) ∂us(k)= ⎣ ⎢ ⎢ ⎢ ⎢ ⎡ ∂y(k) ∂u(k)∂y(k) ∂u(k+1)⋯∂y(k) ∂u(k+s−1) ∂y(k+1) ∂u(k)∂y(k+1) ∂u(k+1)⋯∂y(k+1) ∂u(k+s−1) ⋮ ⋮ ⋱ ⋮ ∂y(k+s−1) ∂u(k)∂y(k+s−1) ∂u(k+1)⋯∂y(k+s−1) ∂u(k+s−1) ⎦ ⎥ ⎥ ⎥ ⎥ ⎤ = T (5.3) Resolving for us(k) gives: us(k)=−[(TT·T+λ· I)−1· TT]·[B · up(k−p)−A · yp(k−p)−ys(k)] (5.4) Where I is an identity matrix. Only are needed the first r (nº of inputs) values of the control sequence ([(TT·T+λ· I)−1· TT]) so the others are unnecessary. Another issue is that to carry out the equation, the desired outputs ys(k) must be given. In this thesis it is taken as 0 and is added an offset in the real output. By this way, taking the first r rows with zero desired outputs leads to the final equation: Us(k)=the first 𝐸 rows of [−(TT·T+λ· I)−1· TT]· ·�B · up(k−p)−A · yp(k−p)� (5.5) 5.2 GPC Parameters Theoretical Effect Before applying the GPC in a system, it is important to understand the effects of its parameters. Basically there are three parameters: the order of the model (p), the weight (λ), and the control horizon (hc), that is related with the prediction horizon. To find the better combination is not trivial, and the better option is to “play” a little with them. The order of the model has been extensively explained on previous chapters. About the control, p can reduce the energy needed to arrive at the desired output and the time to get it. Furthermore increasing this parameter increases very little the computer’s processing time. The weight is a control penalty scalar that limits the control effort and stabilizes the closed-loop system. When λ= 0, the system will be unstable for the most of the structures because the matrix T is rank deficient. On the other hand, a big weight takes more time to get the desired output. Furthermore, the smaller is Predictive Control Algorithm 37 the weight, the bigger is the energy that the input needs. In conclusion, it is needed to select a weight that is a good compromise between the time response of the output and energy consumption. Finally, the growing of control horizon increases the energy required but reduces the time to get the target output. As a very bad point of this parameter, their size is highly related with the time that the computer needs to make the calculations. This is because it works in a similar way of the prediction horizon. The value of hc corresponds to the number of steps that needs the GPC algorithm to drive the output to 0. If for example the control horizon is 1, the input that will get the GPC algorithm would bring the output to 0 only with that input. The problem is that next of that calculated input there is another unknown input. To sum up, a table is presented in the following that shows how some properties of the control algorithm change when its parameters are increased: Table 5.1 Change in the control depending on its parameters. Parameter Energy Needed Time Needed Process time If p − − + If hc + − + + If w − + (no change) 5.3 Energy Consumption Calculus Since the principal objective of this thesis is optimizing the Thermal Vacuum Chamber from an energetic point of view, it is very important to know the amount of energy that will take some system to stabilize itself. The concept of “energy needed” of last section is what is going to be estimated here. The graphics from MATLAB show inputs that seem a curve, but in practise they are not. Actually the inputs are constants in each step. To know the energy consumption is easy, only is needed to do the integral of the input (force) for the differential time. In fact, for the reason told in last paragraph, it is even easier: it is sufficient to compute the area of every input step. These calculi are explained in the next figure: 38 Thermal Vacuum Chamber Control Optimization Figure 5.1 Energy Calculation In the examples that will be shown later and in the application on the Chamber the energy computation results easier yet. Just how the steps have a 1 second duration and the inputs are in the international system units (are used Newtons and Volts), simply is made a sum of every input. 𝐸𝐸𝐸𝐸𝐸𝐸=�𝑢·𝑑𝑑= =𝑢3·(𝑑2−𝑑1)+𝑢2·(𝑑3−𝑑2)+ +𝑢1·(𝑑4−𝑑3) Predictive Control Algorithm Application 39 6. CONTROL ALGORITHM APPLICATION By the same way that is performed with the output prediction, now are shown the tests that have been made with the software about the control. With the same examples that were in the prediction, now is applied the control and is studied how their output and error and energy consumption change with the control parameters. The objective is to prove that the software works well by comparing the results with available data. The developed code for this algorithm can be found in the annex F: Generalised Predictive Control Algorithm. 6.1 Control Example with single Input and Output In the exact code that was used in the section 4.1.2 (Prediction Example with single Input and Output) now has been added and changed code to implement its control. To understand better how works the control, next figure is a case of this example controlled: Figure 6.1 Example 1 controlled with p= hc=2 and w=10. The desired output is 1, and the control applied starts at step 3. The needed control is very low but it can be seen that is not necessary more. 40 Thermal Vacuum Chamber Control Optimization The two next pictures show the mean error and the energy in function of the weight. In addition, the results are obtained with an order model of 3 (blue) and 5 (red). The first figure is about the mean error that has the output respect with the weight applied at the control. The mean error means the average difference between the desired output and the real output along all the process. It is the way used to express the time needed to get the desired output. The bigger is the mean error the bigger is the time that needs the system to stabilizes itself. Figure 6.2 Error in function of weight The first blue line corresponds to an hc=2 the next is with hc=3, until the last blue, whose value is 5. For the red ones is the same, the first is with a control horizon of 2 and the last of 5. This figure accomplishes the theoretic explanation. If the weight is bigger, the system has more attention to the control minimization and less to the error, so it needs more time to get the desired output and the mean error is bigger. In addition, the order model and the control horizon also reduce the error. The following figure shows the energy consumption of the system to get the desired output in function of the GPC weight. Another time the control horizon goes from 5 to 3 steps. Here also works the theory: a big order model and weight consume less energy than littler ones. p=3 hc=2 p=5 hc=5 Predictive Control Algorithm Application 41 Figure 6.3 Energy in function of weight Note that when the weight is little is not clear that the parameters work in the same manner. This is because the system becomes more instable. The weight range use to make the analysis has been from 10-2 until 102. 6.2 Control Example with Multiple Inputs and Outputs Now is performed the example with two inputs and two outputs that has been used in the prediction tests. This time is shown the direct effects that have the control parameters on the input and output of the system. The desired output to get by the control is 1 and 2. Using an order system of 2, a prediction horizon of 2 and a weight of 10-2 are obtained the following inputs and outputs (where YR1 and YR2 are the outputs without control and YP1 and YP2 are the controlled outputs): p=3 hc=5 p=5 hc=2 Bibliography 49 BIBLIOGRAPHY [1] Juang, Jer-Nan and Phan, Minh Q.; “System Identification”, Chapter 10 in Identification and Control of Mechanical Systems, pg. 256-294; Cambridge University Press, 2001. [2] Juang, Jer-Nan and Phan, Minh Q.; “Predictive Control”, Chapter 11 in Identification and Control of Mechanical Systems, pg. 295-331; Cambridge University Press, 2001. [3] Lattuda M., Pavarotti A., Scarpellini N., Trovato D.; Automazione del controllo termico di una camera termovuoto per prove di qualifica spaziale; University tutor: Dr. Bernelli F.; Project manager: Ing. Malnati F.; Facoltà di Ingegneria Aerospaziale, Politecnico di Milano, 2008. ANNEXES A. Pictures about the Chamber Detailed pictures about some components of the Thermal Vacuum Chamber are reported in this Annex. The description of the photos follows: - 1A: Frontal of all the equipment related with the Chamber. - 1B: The Chamber’s back. - 2A: Outside connector gate with the testing object. - 2B: Inside connector gate with the testing object and infrared lamp. - 3A: Infrared lamp detail. - 3B: Two types of thermocouples (centre and top right) - 4A: One of the heating resistances that are in the aluminium plate. - 4B: Pistons pump - 5AB: Chamber’s baseplate design scheme. A B 1 2 3 4 5 B. Alpha and Beta Calculation function [alpha, beta] = arx_batch (u, y, p) % Identification with ARX model in batch mode. % % function [alpha, beta] = arx_batch(u, y, p) % % Function arx_batch identifies an ARX model of order 'p' from arbitrary % input and output data. % % Input parameters: % % u = input histories [l x r], l is the length of the input data % y = output histories [l x m], l is the length of the output data % p = ARX order % % Output parameters: % % The identified model matrix coefficients 'alpha' and 'beta': % % y(k) = - alpha(1) y(k-1) - alpha(2) y(k-2) - ... - alpha(p) y(kp) + % beta(0) u(k) + beta(1) u(k-1) + ... + beta(p) u(k-p) % % alpha = auto regressive coefficients [m x pm] % [ alpha(1) | alpha(2) | ... | alpha(p) ] % % beta = exogenous part coefficients [m x (p+1)r] % [ beta(0) | beta(1) | ... | beta(p) ] % % References: % % [1] Raymond G. Kvaternik, Jer-Nan Juang, Richard L. Bennett % "Exploratory Studies in Generalized Predictive Control for Active % Aeroelastic Control of Tiltrotor Aircraft" % NASA/TM-2000-210552 % % [2] Jer-Nan Juang, Minh Q. Phan % "Identification and Control of Mechanical Systems" % Cambridge University Press, 2001 % number of samples, dimension of output [l m] = size(y) ; % dimension of input r = size(u,2) ; % number of rows of V nrow = r + (m + r) * p ; % allocation V = zeros(nrow, l - p) ; % build first rows of V with inputs V(1 : r, 1 : l - p) = u(p+1 : l, 1 : r)'; % build v auxiliar matrix v = [y'; u'] ; for i = 1 : p ri = r + (r + m) * (i - 1) + 1 ; re = ri + r + m - 1 ; V(ri : re , :) = v( : , p+1-i : l-i ) ; end clear ri re % observer Markov parameters Y = y(p+1:l,:)' * pinv(V) ; % allocate coefficients alpha = zeros(m, m * p) ; beta = zeros(m, r * (p + 1)) ; % beta(0) beta(:, 1 : r) = Y( :, 1 : r ) ; % extract coefficients for i = 1 : p pi = r + 1 + (r + m) * (i - 1) ; pe = pi + (m - 1) ; ci = (i - 1) * m + 1 ; ce = ci + (m - 1) ; alpha(:, ci : ce ) = -Y(:, pi : pe) ; pi = pe + 1 ; pe = pi + (r - 1) ; ci = r + r * (i - 1) + 1 ; ce = ci + (r - 1) ; beta(:, ci : ce) = Y(:, pi : pe) ; end end C. Multistep Matrix Generator with Single input/output function [T, A, B] = matrixs_T_A_B (alpha, beta, s) % This code is the function to obtain the three matrixs (T, A, B) for % multistep prediction. This must be applied only to systems with one % input and one output. % arx order (alpha is m x pm) p = size(alpha,2); %%%% T matrix %%%% T = zeros (s,s); for q = 1:s T (q,q) = beta (1,1); % makes the diagonal end for k = 1:s-1 if k <= p beta_k = beta (1,k+1) ; alphabeta_result = 0; for i = 1:k alpha_i = alpha(1, i); betaki = T (k-i+1,1); % betaki = beta0^(k-i) alphabeta = alpha_i * betaki; alphabeta_result = alphabeta_result + alphabeta; % summation ? end betak = beta_k - alphabeta_result; for q = 1:s-k T (k+q,q) = betak ; end elseif k>p alphabeta_result = 0; for i = 1:p alpha_i = alpha(1, i); betaki = T (k-i+1,1); % betaki = beta0^(k-i) alphabeta = alpha_i * betaki; alphabeta_result = alphabeta_result + alphabeta; % summation ? end F. Generalised Predictive Control Algorithm function [y_pred, u_pred, energy] = pred_cont (u, y, alpha, beta, hc, weight) s =1; %prediction horizon after the control % number of samples, dimension of input [l r] = size(u); %entrada de dades en vector vertical! % dimension of output m = size(y,2) ; %nº columnas % arx order (alpha is m x pm) p = size(alpha,2)/m ; % computation of T A B matrixs: [Tp, Ap, Bp] = matrixs_T_A_B (alpha, beta, s); [T, A, B] = matrixs_T_A_B (alpha, beta, hc); % computation of the other matrixs y_pred = y(1:p,:); %puts the past outputs u_pred = u(1:p,:); for k=p+1 : l-s+1 Yp = y_pred(k-p : k-1,:); Yp = reshape(Yp',size(Yp,1)*size(Yp,2),1); Up = u_pred(k-p : k-1,:); Up = reshape(Up',size(Up,1)*size(Up,2),1); % Calculus of the input control.Formula applied: 11.8 (pag297). %Identity matrix: I = eye(size (T,2)); %Control calculation T_tra = T'; r_row = (-inv((T_tra*T+weight*I))*T_tra); alpha_c = -r_row*A; alpha_c = alpha_c(1:r,:); beta_c = r_row*B; beta_c = beta_c(1:r,:); % cont = -r_row*A*Yp+r_row*B*Up; cont = alpha_c*Yp + beta_c*Up; cont = cont'; %maximum i minimum of the Chamber: (0;3) if max(cont)>3 cont=3; end if min(cont)<0 cont=0; end u_pred (k,:) = cont; Us = u_pred(k : k+s-1,:); Us = reshape(Us',size(Us,1)*size(Us,2),1); % y_pred computation % Ys(k) = T*Us(k) + B*Up(k-p) - A*Yp(k-p) Ys = Tp*Us + Bp*Up - Ap*Yp; Ys = reshape(Ys,m,s)'; %Reestructuration of the matrix %Construction of the complet new output matrix y_pred (k,:) = Ys; end %%% Energy needed %%% %every value of "u" is for 1 second duration, so the integral of the input %in function of the time is only te sum of all input values. u_temp = abs(u_pred); energy = 0; for qq = 1:l energy = energy + u_temp(qq); end end