Modeling and Temperature Control of an Industrial Furnace

(1)

Master of Science Thesis in Electrical Engineering

Department of Electrical Engineering, Linköping University, 2016

Modeling and temperature

control of an industrial

furnace

Hampus Carlborg

Henrik Iredahl

(2)

Master of Science Thesis in Electrical Engineering

Modeling and temperature control of an industrial furnace

Hampus Carlborg Henrik Iredahl LiTH-ISY-EX–16/4978–SE

Supervisor: Urban Johansson

Sandvik AB Fredrik Sandberg

Sandvik AB Kamiar Radnosrati

isy_{, Linköpings universitet}

Examiner: Johan Löfberg

isy, Linköpings universitet

Division of Automatic Control Department of Electrical Engineering

(3)

Abstract

A linear model of an annealing furnace is developed using a black-box system identification approach, and used when testing three different control strategies to improve temperature control. The purpose of the investigation was to see if it was possible to improve the temperature control while at the same time decrease the switching frequency of the burners. This will lead to a more efficient process as well as less maintenance, which has both economic and environmental bene-fits.

The estimated model has been used to simulate the furnace with both the ex-isting controller and possible new controllers such as a split range controller and a model predictive controller (MPC). A split range controller is a control strategy which can be used when more than one control signal affect the output signal, and the control signals have different range. The main advantage with MPC is that it can take limitations and constraints into account for the controlled pro-cess, and with the use of integer programming, explicitly account for the discrete switching behavior of the burners.

In simulation both new controllers succeed in decreasing the switching and the MPC also improved the temperature control. This suggest that the control of the furnace can be improved by implementing one of the evaluated controllers.

(4)

(5)

Acknowledgments

First of all we want to thank Sandvik for the opportunity for us to do this thesis with you. More specifically, we want to thank Urban Johansson and Fredrik Sand-berg at Sandvik for the support during this time.

We want to thank our supervisor Kamiar for the comments on the report. A tremendous thanks to our examiner Johan Löfberg for the advices which made this thesis and report a lot better.

Finally, we want to thank our families, fellows and friends for these five years at Linköpings University.

Linköping, May 2016 Hampus Carlborg and Henrik Iredahl

(6)

(7)

Notation

Abbreviations

Notation Meaning

mpc _{Model Predictive Control}

pid Proportional, Integral, Sifferential (Controller)

dmpc Distributed Model Predictive Control

pem Prediction Error Method

miso Multi Input Single Output

mimo Multi Input Multi Output

(10)

(11)

1

Introduction

The purpose with this thesis was to evaluate the control of the temperature in an annealing furnace, both the existing control and possibly new control strategies. Two main goals were produced, the first one was to reduce how often the burners switch between on and off. The life span of the burner can be increased if the frequency of turning on and off the burners can be decreased. In addition, main-tenance costs will most likely be reduced as well. The second goal was to achieve a smoother temperature of the furnace.

To accomplish this the work was divided into two parts, model estimation and control design. The model estimation part entailed development of an accurate and reliable model of the temperature behaviour in the furnace. The model will be used to compare different control designs and also in the implementation of a model based control design.

1.1 Background

The work for this master thesis has been performed at Sandvik Materials Technol-ogy, in Sandviken. Sandvik is a global industrial company and Sandvik Materials Technology is the business area which mainly develop and manufacture stainless steel products. In the production of steel, furnaces are used for heat treatment. One kind of heat treatment is annealing, where steel bars are heated up and then cooled down fast. The purpose of the process is to get the steel more ductile and reduce the hardness. These properties of the steel are decided by the temperate profile during the annealing process. Small changes of the temperature can in the worst case cause an entire batch of product go to waste, which cost a lot of money. If this happens to many times it could lead to the buyers start losing faith and start searching for other opportunities. Since the steel industry having

(12)

2 1 Introduction tough times nowadays is it more important than ever to have cost effective pro-cesses with good profitability.

The problem today is the temperature distribution in the furnace. For some sec-tions in the furnace the temperature is too high even when the burners are off. Every year the goal is to have maintenance for the whole furnace at a specific week and nothing during the other weeks. Some of the parts that wear the most in the furnace are the burners, as every switch between on and off causes wear on them. By reducing the switching, the maintenance and replacement of burners can hopefully be reduced. The big saving is in that the furnace can operate longer without interruption since it takes a long time for the furnace to cool down and heat up for the maintenance.

1.2 Related work

There exist several reports about modelling an annealing furnace, but there ap-pear to be no standard, since there are differences between the furnaces in the layout, the numbers and kinds of burners, kind and shape of the material, and fuel. Those difference make it difficult to use the same approach from other work of modeling furnaces since the properties of the furnaces differ and the assump-tions are not valid. In [8], a comprehensive mathematical model of an annealing furnace is developed, the model takes both radiation and convective heat transfer in consideration for the components in the furnace. The different components in model are the flue gas, the insulation and the product (strip). In [7], a 3D finite el-ement model is developed using COMSOL Multiphysics software to calculate the 3D temperature distribution from radiative heat transfer. The result are used to improve a simplified model in 2D. In this annealing furnace is the material also a strip. In [1], the developed mathematical model is based radiation heat transfer in the furnace. In the furnace is slab reheated and the temperature field inside the slabs is part of the model. The flue gas in this furnace is non-participating in this model unlike the flue gas in [8]. The main difference compared to the the work presented in this thesis is our use of black-box modelling, compare to the mentioned references where grey/white-box modelling was used.

1.3 Limitations

The data collection was limited since the furnace was used in production. There-fore, it was not possible to specify the input during the data collection and it was done during feedback.

It was not possible during the thesis to either add or modify the position of the burners. No investigation of modifying the burners has been done either.

(13)

1.4 Method 3

1.4 Method

The first part was to do a pilot study and find related work. This part was per-formed in order to find opportunities for both the modeling and the controlling part.

This was followed by data collection in the programIBA Analyser which is used

to log the different signals in the furnace. The collected data was imported into MATLAB to estimate and validate a model of the furnace.

The last part was to take the model into SIMULINK to test different controllers. Three different control strategies were applied, Split range, model predictive con-trol (MPC) and distributed model predictive concon-trol (DMPC). The MPC prob-lems were formulated by using the toolbox YALMIP [5] and were solved with the solver MOSEK [2]. All strategies were then evaluated, in terms of temperature and numbers of switches for the burners.

1.5 Thesis outline

The thesis is organized as follows: • Chapter 2 presents the furnace.

• Chapter 3 gives the background to the control design. • Chapter 4 presents the model estimation.

• Chapter 5 presents the control design and the result of the simulation. • Chapter 6 presents what can be done in the future based on the results from

(14)

(15)

2

System Description

The furnace that has been analysed during this thesis is an annealing furnace. Annealing is a process in which the material (steel bars in our case) heats up to a certain temperature to obtain specific material properties and then cools down fast to capture the properties. This type of process is well known and common for steel products. Since the temperature profile determines which properties the material will have, it is important that the furnace achieves a given reference tem-perature. A too big temperature difference can cause a whole batch of processed material to be wasted. The furnace was divided into three different temperature zones, where the reference temperature is highest for the zone where the bars exits and lowest where the bars enter the furnace.

The work order specifies how long the bars should be in the furnace and the reference temperature for the batch. A walking beam will push the bar through the temperature zones. In connection with the walking beam is used, the doors opened before / after to take out / in bars from the oven. The material that goes into the furnace were cylindrical steel bars so the material could rotate in the fur-nace.

For this thesis the furnace for heating up the steel bars has been analysed and not the cooling process. Since the furnace has been in use during the thesis, all the data has been collected under real production. This made it impossible to perform any specific tests.

(16)

6 2 System Description

2.1 Furnace

The furnace is divided into three different temperature zones, where the first zones had the lowest reference temperature and the last zone had the highest reference temperature. To achieve a larger temperature difference between zone 1 and 2 a separator was mounted in the roof between the zones. Every temperature zone was divided into three sections. It was necessary to divide the zones into sections since the heat flow between sections in the same zone was big and cause temperature differences inside the zones. If the temperature had been the same for each zone the sections would not have been necessary. Figure 2.1 shows a sketch of the furnace. The arrows indicate where the bars enter and exit the furnace. It can be seen that the bars enter into the furnace through zone 1 west and exit from zone 3 west. Every time a bar enter or exit, a door needs to opened which causes a temperature difference for some sections. When the bars enter the furnace it has the same temperature as the room.

Figure 2.1:This is a sketch of the layout of the furnace seen from above, the

arrows indicate where the bars enter and exit, the flames indicates where the burners are placed. The numbers indicate the number of the sections.

The width of the furnace is about 7 meters, length 9.25 m and height 3.3 m. The bars are transported horizontally into the west section of zone one, and depending on the length of the bars it could occupy one to all three sections. All the burners for each zone is in line with each other. It would be a too big effort to change the place for the burners so the numbers and the position of them can be seen as constant. The furnace has 27 compartments, which each can contain one bar. For bars with large diameters they cam only be placed in every second compartment. The bars are transported to the next compartments by a walking beam. Figure 2.2 shows the measured temperature for a batch. It is easy to see the zig-zag behaviour in the temperature, with a period of roughly 750 seconds, caused when the walking beam are transporting the bars forward.

(17)

2.1 Furnace 7 0 500 1000 1500 2000 2500 3000 3500 940 960 980 1000 1020 1040 1060 1080 1100 1120 Original Temperatures Temperature (celsius) Time (seconds) Z1 W Z1 M Z1 E Z2 W Z2 M Z2 E Z3 W Z3 M Z3 E

Figure 2.2:The measured temperature for a certain order. Note the zig-zag

behaviour in the temperature, with a period of roughly 750 seconds, caused when the walking beam are transporting the bars forward.

In sections west and east in zone 1 there is an exhaust for the flue gas which caused heat transfer from zone 3 into zone 1. The heat streams have a large im-pact on the temperature in zone 1 which cause the temperature to remain too high even when the burners are off in east and west zone 1. For each section there was different numbers of burners, varying from one to four as illustrated in Figure 2.1.

Every section of the furnace had a PI-controller that controls the temperature in the section. The burners had two modes, being either on or off. The control signal from the PI-controller gave a percentage on which capacity the burners in the section would operate on. The sensor for measuring the temperature in every section was mounted 0.4 meters from the roof. When a new order with a different temperature was requested it could take a couple of hours to set it up depending of how large the temperature difference will be. The temperature ref-erence is ramped slowly to the new level. The main point of the slow change of the reference is to ensure that the whole section has the right temperature and not only close to the sensor.

(18)

8 2 System Description

2.2 Burner system

The burner system had a separate control system consisting of PI-controllers. The number of burners vary between sections, e.g. in zone 1 Mid the PI-control con-trols four burners and in zone 2 West it was only one burner. The controlled burners had two modes either on or off. Since the burners only operated with two modes the control signal was converted to a burning time for the burners. The PI-controllers were implemented as series PI-controllers.

F(s) = K(1 + I ∗1 s)

The parameters were K = 3.36 and I = _93.461 . The burning time was given by

the cycle time multiplied by the control signal and the numbers of burners in the section. To minimize the frequency of switching between on and off the burning time was set to zero if it was less than 5 seconds and to maximum burning time if it was greater than maximum burning time minus 10 seconds. The burner system can be summarized in the following algorithm:

Algorithm for burner system:

1. Wait until a cycle time have passed.

2. Calculate the burn time as the cycle time multiplied withthe control signal and the number of burners in the section.

3. If the burn time is less than 5 seconds set it to zero.

4. If the burn time is greater than maximum burn time minus 10 sec-onds set it to maximum burn time.

5. Set the next burner to burn for the burn time.

6. Update next burner, whose burning time shall be calculated the next time, to the last burner to have it burn time set.

7. Repeat from step one.

If a section had two burners and the control signal was 50 %, the burners burned on each half of the time, or if there were two burners with the signal 25 % half of the time no burner would burn. Figure 2.3 shows the burner system for section 2, which has 4 burners, with the control signal at 25 %, the burners operates in series, when the first one is done the next starts and so on. Figure 2.4 shows the same section for 60%, the burner are now working more in parallel. In the system today none of the burners in zone 1 and 2 east are on since it is already to hot there. There is small delay before the burner starts to burn.

(19)

2.2 Burner system 9 0 50 100 150 200 250 300 350 Time -0.50 0.51 1.5 Burner 1 0 50 100 150 200 250 300 350 Time -0.50 0.51 1.5 Burner 2 0 50 100 150 200 250 300 350 Time -0.50 0.51 1.5 Burner 3 0 50 100 150 200 250 300 350 Time -0.50 0.51 1.5 Burner 4

Figure 2.3:How the burner system act with an input of 25%, a one represent

(20)

10 2 System Description 0 50 100 150 200 250 300 350 Time -0.50 0.51 1.5 Burner 1 0 50 100 150 200 250 300 350 Time -0.50 0.51 1.5 Burner 2 0 50 100 150 200 250 300 350 Time -0.50 0.51 1.5 Burner 3 0 50 100 150 200 250 300 350 Time -0.50 0.51 1.5 Burner 4

Figure 2.4:How the burner system act with an input of 60%, one represent

(21)

3

Control design theory

The purpose of this chapter is to give the background for the different control designs that were implemented in this thesis. The different controllers that were used were split range control, MPC and DMPC.

3.1 Split range control

Split range control is a simple control strategy which can be used when there are more than one input that controls the same output. The basic idea is to divide the control signal into different ranges for the input signals. In the case of two burners, only one will burn in the range 0-50 % and the other is off, and in the range 50-100 % one will burn all the time and the other is burning according to the control signal. A setup for a system with a split range controller can be seen in Figure 3.1 Burner 1 Burner 2 PI-Control Split-range ref temp % u1 u2

Figure 3.1: A sketch for the setup of the system with the split range

con-troller implemented.

A split range control can decrease frequency the burners turn on and off since 11

(22)

12 3 Control design theory

all the burners in a section, except one is burning at 100 or 0 %. The split range control have two parameters, swap time and cycle time. Swap time decides the time between the swap of the range for the burners. The burners have different positions inside the section, and therefore affect the temperature in its own sec-tion and nearby secsec-tions differently.

The use of swap time is to minimize the effect of the burners position such as uneven temperature within the section. The cycle time is used in the same way as in the current burner system, i.e., to convert the control signal to a burn time by multiplying the control signal with the cycle time. An example of the signal from the burners can be seen in Figure 3.2. The split range controller can be summarized with the following algorithm:

Algorithm for split range:

1. Wait until it is time to calculate the burn time.

2. If the time that has passed since the last swap is greater than the swap time.

• Change the operating range for the burners. 3. Calculate and set the burn time for each burner. 4. Repeat from step one.

(23)

3.2 Model Predictive control 13 0 50 100 150 200 250 300 350 −0.5 0 0.5 1 1.5 Burner 1 Time 0 50 100 150 200 250 300 350 −0.5 0 0.5 1 1.5 Burner 2 Time 0 50 100 150 200 250 300 350 −0.5 0 0.5 1 1.5 Burner 3 Time 0 50 100 150 200 250 300 350 −0.5 0 0.5 1 1.5 Burner 4 Time

Figure 3.2:This shows how the burners act when the split range controller

have an input signal of 60 %, 1 represent that the burner is on and 0 that is off.

3.2 Model Predictive control

In this section the MPC is introduced. We first present the background and a ba-sic MPC formulation. An advantage of the MPC is its flexibility and that the baba-sic formulation can be extended to include reference tracking, feedforward control and relaxed constraints, as described. For further reading see [6].

3.2.1 Background

Model predictive control is based on solving an optimization problem for each sample instant. The solution to the optimization problem is a sequence of control signals to be used now and in the future. Of the sequence of control signals, only the first is applied to the controlled system and the rest is discarded. This procedure is repeated for each sample instant and achieves a feedback control as the optimization problem is performed when new measurements are made available. If the whole sequence would be applied it would be open loop control. The model of the controlled system is a part of the optimization problem and is used to predict the future states. The MPC can be seen as it performing an open loop optimization for each sample instant.

(24)

14 3 Control design theory

3.2.2 Description of the model predictive controller

This section will introduce the model predictive controller with a theoretical ex-ample. The objective for the MPC is to bring the system states and the input to zero, subject to the system dynamics and limitations. In the basic setup the cost function is quadratic and the constraints are linear.

minimize u(·) N X j=1 ||_{x(k + j)||}2 Q1+ ||u(k + j − 1)|| 2 Q2 subject to x(k + i) = Ax(k + i) + Bu(k + i) xmin≤xi ≤xmax umin≤ui ≤umax (3.1)

Here N represents the prediction horizon which determine how many steps

for-ward the controller predicts. The norm ||x||2_Qis the euclidean norm weighted by

the weight matrix Q, xTQx. The difference between the weight matrices Q1and

Q2 determinate the behaviour for the controller. If Q1is large in comparison to

Q2 will the system states will converge to zero fast, at the cost of large control

inputs. If Q1is small in comparison to Q2will the control inputs be small, at the

cost of that, the system states will converge to zero slow.

The equality constraints represent the system dynamics from the model of the controlled system. The inequality constraints are bounds for the system states and/or the control signal, which are typically for safety regulations or saturation. In summary the MPC can be explained by the following algorithm:

Algorithm for MPC:

1. Measure the system state x(k)

2. Compute the control signal sequence u( · ) for the problem in (3.1) 3. Use the first element u(k) in the sequence from the previous step

during one sample 4. Time update, k = k + 1 5. Repeat from step 1

(25)

3.2 Model Predictive control 15

3.2.3 Reference Tracking

In the formulation in (3.1) the system will be driven to a steady state at the origin. In order to steer the system to any other state than the origin, the MPC formula-tion must be extended to include reference tracking. This can partially be accom-plished by minimizing the difference between a state and reference. This gives a a conflict in the MPC formulation which is to have the states follow a reference and try to keep the input zero. One solution is to minimize the increment of the input signals instead of its amplitude. The reformulation (3.2) gives the MPC controller reference tracking.

minimize u(·) N X j=1 ||_{x(k + j) − r(k + j)||}2 Q1+ ||u(k + j − 1) − u(k + j − 2)|| 2 Q2 subject to x(k + i) = Ax(k + i) + Bu(k + i) xmin≤xi ≤xmax umin≤ui ≤umax (3.2)

3.2.4 Feedforward Control

In feedforward control the disturbance is measured and with used in the control computation to remove some of the disturbance effect on the system. It thus has an advantage compared to feedback control which has to wait for the disturbance to effect the system before it can take suitable action. The effectiveness of feedfor-ward control is limited by the model of the measurable disturbance to the output, and any remaining unmeasurable disturbances. Therefore feedforward control is implemented together with feedback control. The feedforward control part re-moves most of the measurable disturbance effect and feedback control take care of the rest and also any unmeasured disturbance. Feedforward control is imple-mented in MPC by including the disturbance effect in the predicted future out-puts, which can be seen in (3.3). The optimizer then accounts for the disturbance when it computes the control signal.

(26)

16 3 Control design theory minimize u(·) N X j=1 ||_{x(k + j)||}2 Q1+ ||u(k + j − 1)|| 2 Q2 subject to x(k + i) = Ax(k + i) + Bu(k + i) + Bddm(k + i) xmin≤xi ≤xmax umin≤ui ≤umax (3.3)

3.2.5 Relaxed Constraints

A serious problem for the MPC is when the optimization problem is infeasible. This can happen because an unexpected large disturbance occurred and it is im-possible for the controller to prevent that a state breaks its constraint. Therefore it is important that the MPC controller has a strategy to handle infeasible prob-lem. A systematic strategy is to make relaxed constraints which can be broken if there is no admissible solution [6]. The constraints on the input signal are usually "hard" and cannot be broken since it commonly represent a physical limitation on the actuator. Therefore it is usually the constraints on the output which are relaxed.

A constraint is relaxed by introducing a new variable which is non-zero when the constraint is violated. Such a variable is called slack variable. To avoid viola-tion of constraints when unnecessary, a penalty term is added to the objective in order to try to obtain a zero slack. The optimization problem then becomes

minimize u(·) N X j=1 ||_{x(k + j)||}2 Q1+ ||u(k + j − 1)|| 2 Q2+ λ|| subject to x(k + i) = Ax(k + i) + Bu(k + i) xmin− ≤ xi ≤xmax+  umin ≤ui ≤umax 0 ≤  (3.4)

The reason for not using the euclidean norm on in the cost function is that all active constraint will be violated for all finite λ, even when it is not necessary.

From the true constrained x∗there exist a move to an infeasible solution x∗+ d

(27)

3.3 DMPC 17

norm will the cost of violating the constraint be O(2) which gives a net

reduc-tion in the cost for small . With the 1-norm in the cost funcreduc-tion a sufficient λ will ensure that violation occur only when the original problem is infeasible.

3.2.6 Complete model

The different modifications of the basic formulation can be used together and gives the following optimization problem, which will be the form used later.

minimize u(·) N X j=1 ||_{x(k + j) − r(k + j)||}2 Q1+ ||u(k + j − 1) − u(k + j − 2)|| 2 Q2+ λ|| subject to x(k + i) = Ax(k + i) + Bu(k + i) + Bddm(k + i) xmin− ≤ xi ≤xmax+  umin≤ui ≤umax 0 ≤  (3.5)

3.3 DMPC

Another way to use the function of MPC is to enlarge the system with more MPC controllers. This type of control strategy is called distributed model predictive control. The meaning of DMPC is to split the problem into sub systems, where every MPC controller computes a control input locally. This will make it easier for the solver, as the optimization problems will be smaller. However, when do-ing this, the interaction between the different controllers and parts of the systems will not be part of the models used locally in every MPC controller. One approx-imation, which we will use, is to assume that the states from other zones will be constant over the prediction horizon. For further reading about DMPC, see [9].

(28)

(29)

4

Model estimation

This chapter is about the identification and validation of the model for the system described in Chapter 2. The workflow for the model estimation was an iterative process; it started with a simple model that was extended until no improvements could be found. The main idea with the model was to keep it simple with rele-vant relations, which capture the dominant behaviour of the system.

4.1 Black-box modelling

There are three different approaches of deriving models; White-box, Grey-box and Black-box. The White-box model is based only on known information about the system. If all information about the system is known and the model has the right structure, the model and the system will be identical. In the case of the fur-nace much of the information is missing, such as the heat transfer between the sections and heat losses through the walls. It could be possible to model those relations based on thermodynamics, but many parameters are unknown such as the flow of the flue gas and need to be estimated. A White-box model of the fur-nace would be very complex and some relations are hard to describe with physics. For example the temperatures zig-zag behaviour caused by the walking beam. The advantage with a Black-box model is that it can be used when no informa-tion about the system is known. The Black-box model describes the relainforma-tions between the input and output signals without having a structure based on phys-ical relations of the system. The model can be in the form of equations, tables or graphs. A possible problem could be that the model only performs for precisely those cases for which it was derived.

(30)

20 4 Model estimation

The Grey-box model can be seen as a hybrid between Black-box and White-box model. The part of the system whose information is known is modelled as a White-box and the rest of the system as a Black-box. This model was neglected with the same reason as the White-box, it is too hard to model the thermodynam-ics.

The Black-box model was chosen because it is easier to implement and the model is a tool in this thesis. The possible improvements by using a Grey-box was con-sidered low, as said before a large part of the system is unknown.

4.2 Discrete-time state space

A broad definition of a system is an object or a collection of objects whose prop-erties are observed. A model can be used as a tool to examine the system without doing any experiment on it. To experiment with a new controller in the furnace is ineffective and will cost money in the lost production. The advantage with a model is that it gives us the possibility to try different control strategies without implementing them in the furnace.

A common way to represent a system is the model in the following form

˙x = f (x, u)

y = h(x, u) (4.1)

where u is the input signal, y is the output signal and x is called the state. The equations describe the whole system’s behaviour. This structure is called state space. Specializing to the linear case, and discretizing the continuous-time dy-namics using, e.g., a first-order hold, gives us a discrete-time state space model

x(kTs+ Ts) = Ax(kTs) + Bu(kTs)

y(k) = Cx(kTs) + Du(kTs) (4.2)

where the matrices A,B,C and D describes how the future state depends on the current state and the input as well as how the output depends on the state and the input. For the furnace, the states are the temperatures; the input signals are the burners, the position of the walking beam and the position (open/close) of the doors. The output signals are the measurements of the temperatures. Since the furnace is an unknown system the objective is to estimate the unknown pa-rameters in the matrices. The furnace itself is a nonlinear system and therefore the linear model is only a rough approximation of the furnace.

(31)

4.3 Data collection 21

4.3 Data collection

The data of the furnace is collected during production. Thus the data was col-lected during feedback. For the furnace there is a logging system called IBA ANALYSER which logs all sorts of data. The collection lasted for 1.5 hours with a sampling rate of 5 seconds. There exist different methods of handling data col-lection during feedback such as the direct approach, the indirect approach and the joint input-output approach [3]. In this thesis, the direct approach was cho-sen where the effect of the feedback is ignored. Since the priority was the control design and the model as said before was a tool.

A problem with the data collection was that the temperatures were too high for some sections. The burners in those sections did not burn, and thus there is no in-formation in the collected data set on the relation between these burners and the resulting temperature. Those burners are still included in the model but their parameter values are just an approximation. The reason for including approxi-mate values of these parameters anyway was to enable us to use these burners in simulation.

4.4 Parameter estimation

To estimate the parameters for a model one of many methods is the prediction error method [4]. The vector θ includes the parameters that will be estimated,

in our case A,B,C,D . It will be used by the prediction ˆy(t|θ) which depends on

y(t-1). The prediction error (t, θ) is the difference between the output y(t) and ˆ

y(t|θ).

(t, θ) = y(t) − ˆy(t|θ) (4.3)

The loss function for a data set with N samples is written as

VN(θ) = 1 N N X t=1 2(t, θ) (4.4) ˆ

θN is given as the parameter which minimizes the loss function

ˆ

θN = arg min

(32)

22 4 Model estimation

To perform the prediction error method in MATLAB we used PEM from the SYSTEM IDENTIFICATION TOOLBOX. The input to the function is the estima-tion data from the system and the model with the parameters that should be estimated. P EM has an option on estimation focus which weigh the loss function for specified frequencies. Two different focuses were used, prediction and sim-ulation. For prediction has a loss function which minimizes the one step ahead prediction. For simulation which has a loss function that favours the frequency range where the input has the most power.

4.4.1 Simplification and constraints

The state space model in (4.2) is a multiple input multiple output (MIMO) system which made the parameter estimation difficult since it was an optimization prob-lem with many parameters. To make it easier to optimize the system was split into nine different multiple input single output (MISO) systems, then the MISO-systems were converted back to a MIMO-system. The advantage of converting it back was that it was easier to implement in SIMULINK and include it in the MPC formulation. The MISO-systems have only one state, which is the temper-ature in its section. The input signal was expanded to include the tempertemper-atures from nearby sections. An example is given below for how section one was treated as a MISO-system.

T1(k + 1) = aT1+ [b1 · · · b7][T2 T4 B1 B2 W D1 D2]T

T1 is the temperature in section 1, T2 in section 2, T4 in section 4, B1 is the

burner in section 1, B2is the burner in section 2 which is closest to section 1, W

is the walking beam, D1 the door in section 1 and D2the door in section 9. The

parameters for the MISO system were then included into the corresponding place in the A and B-matrix in the MIMO system, which is given in an example below.

A1,1:9= [a b1 0 b2 0 0 0 0 0 0]

B1,1:20 = [b3 b4 0 · · · 0 b5 b6 b7]

To avoid running into a local minimum, constraints were added to the parameters that have a physical interpretation. The estimation was first done without the constraints, but then some parameters were negative for the burners, which had caused a temperature decrease when they were on. A burner which is burning can not decrease the temperature. All the bounds of the parameters can be seen in Table 4.1. The parameters that involve the doors and the walking beam have no constraints since it is hard to interpret its effect on the temperature. When a door is opened there can be an increase in heat transfer between two sections, the temperature in one section will increase and the other will decrease.

(33)

4.5 Model Validation 23

Table 4.1:The constraints for the parameters.

Name Lower bound Upper bound

Temperature from their own section 0.7 ∞

Temperature from burner in section 0.1 ∞

Temperature from burner in nearby section 0 0.4

4.5 Model Validation

To validate how good the models were three different methods were used. The first method was cross validation, which gives a percentage of how well the out-put of the model represents the outout-put from real system. Number two was resid-ual analysis, which shows the correlation between input and prediction error and also the autocorrelation of the prediction error. The third method was the k-step prediction; it shows how good the model predicts the future for the k steps ahead.

4.5.1 Cross validation

Two thirds of the data was used for estimation while the remaining third was used for validation. The use of different data in validation and estimation is called cross-validation. It is used to prevent overfit where the model would de-scribe the noise instead of underlying relations.

One way to validate a model is to compare its output with the measured input against the measured output. The model fit, in percentage, is given below

Mf = (1 − (t,θ)ˆ y (t)−y(t)¯ ) ∗ 100 (4.6)

(t, ˆθ) is the model error from (4.3), y(t) is the measured output and ¯y(t) is its

average, || · || is the L2 norm. The model fit is an easy way to compare different

models. The model fit is a measure of how well the model represents the real system.

4.5.2 Residual analysis

Another way to validate the model is to look at the residuals of the model. One residual is the cross correlation between prediction error and the inputs. In best case the input signal would be independent of the prediction error, if not the model is missing some dynamics from the system.

(34)

24 4 Model estimation ˆ Ru(τ) = 1 N N X t=1 (t + τ)u(t), |_{τ| ≤ M} _(4.7)

The auto correlation of the prediction error is given as,

ˆ R(τ) = 1 N N X t=1 (t)(t + τ), |_{τ| ≤ M} _(4.8)

If the auto correlation is outside the 99% confidence interval for some τ then the output depends on the output τ step before.

4.5.3 k steps prediction

The k-step prediction is useful to test the model predictive power k steps ahead in time. The model output is predicted at k steps ahead based on the information

at time k0, which is given by

ˆ

Xk0+k+1|k0= A ˆXk0+k|k0+ BUk0+k, k = 0, 1, ...

ˆ

Yk0+k|k0= C ˆXk0+k|k0

(4.9)

The predicted output ˆYk0+k|k0is compared against the measured output Yk0+kthe

same way as in the cross validation in (4.6) where is replaced. In this thesis a predictive controller is used and therefore it is important that it can predict the signals and thereby choose the best control signal.

4.6 Final model

In this section the structure of the final model is presented. The A-matrix repre-sents how future state depends on the current state, in this case the temperature. The diagonal in the A-matrix determines how much temperature every section contains from the previous sample, i.e., the decay in temperature given no in-flux or loss of heat to surrounding sections, or addition of heat using burners. The other parameters represent the effect of the temperature from other sections. The final values can be found in Appendix A.1.

(35)

4.6 Final model 25 A =                                   a1,1 0 0 a1,4 0 0 0 0 0 a2,1 a2,2 0 0 a2,5 0 0 0 0 0 a3,2 a3,3 0 0 a3,6 0 0 0 a4,1 0 0 a4,4 a4,5 0 a4,7 0 0 0 a5,2 0 a5,4 a5,5 a5,6 0 a5,8 0 0 a6,2 a6,3 0 a6,5 a6,6 0 0 a6,9 0 0 0 a7,4 0 0 a7,7 a7,8 0 0 0 0 0 a8,5 0 a8,7 a8,8 a8,9 0 0 0 0 0 a9,6 0 a9,8 a9,9                                   (4.10)

The B-matrix represents the effect on the temperature from the burners, the doors and the walking beam. The values are found in Appendix A.2 for the different models. BT =                                                                                   b11 0 0 0 0 0 0 0 0 b21 b22 0 0 0 0 0 0 0 0 b32 0 0 0 0 0 0 0 0 b42 b43 0 0 0 0 0 0 0 b52 b53 0 0 0 0 0 0 0 0 b63 0 0 0 0 0 0 0 0 0 b74 0 0 b77 0 0 0 0 0 b84 b85 0 0 0 0 0 0 0 0 b95 0 0 b98 0 0 0 0 0 b10,6 b10,7 0 b10,8 b10,9 0 0 0 0 0 b11,6 0 0 0 0 0 0 0 0 0 b12,7 0 0 0 0 0 0 0 0 b13,7 b13,8 0 0 0 0 0 0 0 0 b14,8 0 0 0 0 0 0 0 0 b15,8 0 0 0 0 0 0 0 0 b16,8 b16,9 0 0 0 0 0 0 0 0 b17,9 −_b_18,1 −_b_18,2 −_b_18,3 _b_18,4 _b_18,5 ₀ ₀ ₀ ₀ b19,1 −b19,2 0 0 b19,5 0 0 0 0 −_b_20,1 −_b_20,2 ₀ −_b_20,4 −_b_20,5 ₀ ₀ ₀ ₀                                                                                   (4.11)

Since every state could be measured the C-matrix only will apply as a unit ma-trix.

(36)

26 4 Model estimation C =                                   1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1                                   (4.12)

4.6.1 Validation

(37)

4.6 Final model 27 200 400 600 800 1000 1200 1400 1600 1800 960 965 970 975 980 985 990 995 y1

Validation data (y1) sys1_prediction: 63.18% sys1_simulation: 55.77% 12-Step Predicted Response Comparison

Time (seconds)

Amplitude

(a) Validation for zone 1 west. Notice the models’ smaller amplitude.

200 400 600 800 1000 1200 1400 1600 1800 940 945 950 955 960 965 y1

Time (seconds)

Amplitude

(b)Validation for zone 1 middle. Notice the difference between the two models.

200 400 600 800 1000 1200 1400 1600 1800 998 1000 1002 1004 1006 1008 1010 1012 1014 1016 y1

Time (seconds)

Amplitude

(c)Validation for zone 1 east. The models have the peaks at the right time, but not the right amplitude.

(38)

28 4 Model estimation 200 400 600 800 1000 1200 1400 1600 1800 1056 1058 1060 1062 1064 1066 1068 1070 y1

Time (seconds)

Amplitude

(a)Validation for zone 2 west. Notice the mod-els’ good performance.

200 400 600 800 1000 1200 1400 1600 1800 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 y1

Time (seconds)

Amplitude

(b)Validation for zone 2 middle. The models miss some of the peaks.

200 400 600 800 1000 1200 1400 1600 1800 1063 1064 1065 1066 1067 1068 1069 1070 y1

Time (seconds)

Amplitude

(c)Validation for zone 2 east. Notice that the models have the peaks at the right time, but not the right amplitude

(39)

4.6 Final model 29 200 400 600 800 1000 1200 1400 1600 1800 1093 1094 1095 1096 1097 1098 1099 1100 1101 y1

Time (seconds)

Amplitude

(a)Validation for zone 3 west. Notice the small variation of the temperature.

200 400 600 800 1000 1200 1400 1600 1800 1095 1096 1097 1098 1099 1100 1101 y1

Time (seconds)

Amplitude

(b) Validation for zone 3 middle. The models have the peaks, but not at the right amplitude.

200 400 600 800 1000 1200 1400 1600 1800 1094 1095 1096 1097 1098 1099 1100 1101 y1

Time (seconds)

Amplitude

(c)Validation for zone 3 east. Notice the simi-larities between the models.

(40)

30 4 Model estimation 200 400 600 800 1000 1200 1400 1600 1800 960 965 970 975 980 985 990 995 y1

Validation data (y1) sys1_prediction: 61.33% sys1_simulation: 61.05% Simulated Response Comparison

Time (seconds)

Amplitude

(a) Validation for zone 1 west. Notice the models’ lack of the smaller peaks.

200 400 600 800 1000 1200 1400 1600 1800 940 945 950 955 960 965 y1

Time (seconds)

Amplitude

(b)Validation for zone 1 middle. Notice how the models’ performance gets worse over time.

200 400 600 800 1000 1200 1400 1600 1800 998 1000 1002 1004 1006 1008 1010 1012 1014 1016 y1

Time (seconds)

Amplitude

(c)Validation for zone 1 east. Notice that the models have the peaks at the right time, but not the right amplitude

(41)

4.6 Final model 31 200 400 600 800 1000 1200 1400 1600 1800 1056 1058 1060 1062 1064 1066 1068 1070 y1

Time (seconds)

Amplitude

(a)Validation for zone 2 west. Notice the mod-els’ good performance.

200 400 600 800 1000 1200 1400 1600 1800 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 y1

Time (seconds)

Amplitude

(b) Validation for zone 2 middle. The models lack some of the peaks.

200 400 600 800 1000 1200 1400 1600 1800 1063 1064 1065 1066 1067 1068 1069 1070 y1

Time (seconds)

Amplitude

(c)Validation for zone 2 east. The models have the peaks, but not at the right amplitude.

(42)

32 4 Model estimation 200 400 600 800 1000 1200 1400 1600 1800 1093 1094 1095 1096 1097 1098 1099 1100 1101 y1

Time (seconds)

Amplitude

(a)Validation for zone 3 west, notice the small variations of the temperatures.

200 400 600 800 1000 1200 1400 1600 1800 1095 1096 1097 1098 1099 1100 1101 y1

Time (seconds)

Amplitude

(b)Validation for zone 3 middle, both models have the similar model fit.

200 400 600 800 1000 1200 1400 1600 1800 1095 1096 1097 1098 1099 1100 1101 y1

Time (seconds)

Amplitude

(c)Validation for zone 3 east, the models cap-ture the peaks, but not with the same ampli-tude.

(43)

4.7 Discussion 33

4.6.2 Residuals

Below two of the eighteen residuals is presented, since the residuals have almost the same behavior the rest are found i Appendix C.

0 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 e@y1 AutoCorr 0 XCorr (u1) 0 XCorr (u2) 0 XCorr (u3) 0 XCorr (u4) 0 XCorr (u5) 0 XCorr (u6) 0 XCorr (u7) 0 XCorr (u8) 0 XCorr (u9) 0 XCorr (u10) Residue Correlation Lag Amplitude

(a)The residual for Zone 1 Middle with pre-diction model. 0 -0.2 0 0.2 0.4 0.6 0.8 1 e@y1 AutoCorr 0 XCorr (u1) 0 XCorr (u2) 0 XCorr (u3) 0 XCorr (u4) 0 XCorr (u5) 0 XCorr (u6) 0 XCorr (u7) 0 XCorr (u8) 0 XCorr (u9) 0 XCorr (u10) Residue Correlation Lag Amplitude

(b)The residual for Zone 2 Middle with pre-diction model.

Figure 4.7:The figures show two residuals with estimating from prediction.

4.7 Discussion

It is hard to model the furnace with thermodynamics since it will be a very com-plex model and information about the system is missing. Hence, we chose to model the furnace as a Black-box model. One main advantages of a Black-box model is its easy implementation. We believe that the difference between a Black-box and a Grey-Black-box in model quality are small since much information is

(44)

un-34 4 Model estimation

known. We chose to model the system as a State-space linear model since we thought it would capture the most important aspects of the system. For a MIMO system, it can be difficult to determine which inputs that influence which out-puts, especially when the inputs are binary and can be similar. Therefore, to ensure that the model had only relevant relations, a simple model was chosen. A risk with a too advanced model is that it performs better in a specific case, but not in the general case. For the parameter estimation some restrictions were added. The point of this was to ensure reasonable parameters, such as enforcing that burners only can increase temperature

The effect of the temperature when the walking beam is in use is strange and hard to explain. It causes a zig-zag behaviour in zone 1 but also in zone 2 but of smaller degrees. The peaks can be caused by an increase in the heat radiation from below or that the heat flows are changed. We think the drop is caused when the walking beam is in the end position and air is leaking in from below.

A large difference between the model and the furnace is that the furnace is a nonlinear system. The temperature in the model is smoother and it lacks some of the peaks and drops. In general the peaks have smaller amplitude in the model. The heat losses in the model is seen as the diagonal in the A-matrix since A < 1. As mentioned before the temperature was too high in three sections during the data collection. We chose to include the burners for those sections in the model, but we could not estimate the parameters for them. With simulation with other reference temperature or different controller it could be a case when these burn-ers should be on. The parameter values for the burnburn-ers are approximations based on the parameters for other burners.

The model fit for some sections were really low for the different models. Some reasons can be that the furnace is a nonlinear system and that the data collection was done during feedback. The problem with closed loop data is the correlation between the input and the noise, which typically cause bias in the parameter es-timation. For zone 3 it was hard to model the temperatures since it was close to the given reference temperature and the large part of variation in signal could be caused by noise. Therefore, it was more important to focus on good models of the other zones, since those zones have the biggest opportunity to improve the temperature control. The result from the residual analysis were similar for both models with the worst result in the cross correlation with the burners, it suggests that it exists more information about the system from the burners. In the end there was a small difference between the models and the prediction model was chosen since the model was better for zone 1, as the potential for improvement is greatest in the zone first.

To improve the model one opportunity could be to include the temperature of the bars in the model. The bars had probably have a large impact on the temperature in the furnace since they will cause the temperature in the furnace to decrease

(45)

4.7 Discussion 35

when they are heated. The temperature of the material during its time in the fur-nace could perhaps be estimated with a Kalman filter. Today, the temperature of the material is measured at the end of the furnace. Another improvement could be if the model includes the exhausts, as this could give an advantage of how the heat flow is moving from zone 3 to 1. If this information was known, the dif-ference between the Black-box and the Grey-box could be larger, and Grey-box could be the better option.

Some other improvements that will be harder to integrate could be to add more temperature sensors. This will give more information for the estimation. If some sensors are mounted between the zones it could make it easier to estimate a better heat exchange between the zones. With more sensors in each section, it will give a better picture of temperature in it. The sensor configuration used today does probably not give a complete picture of the whole section. Another way could be to perform specific identification experiments on the furnace. For example, have the same temperature for the whole furnace, and try different burners and see their effect. This would more precisely show the effect of each burner on the temperature in the system and thus most likely allow for a better model.

(46)

(47)

5

Control design

This chapter will show the final controllers and which specification that was given for the controllers.The different controllers were simulated with the model that was derived in Chapter 4. The result from the simulation will be presented with both tables and figures. Finally, the chapter ends with a discussion about the different controllers and their result in the simulation.

5.1 Specification

The aim with a new controller design was to decrease the frequency in which the burners switch between on and off as well as a better temperature control. There are no specific requirements on the frequency of the switching, but any reduction would increase the life span for the burners, which allows them to be used longer before replacement. For the temperature control it is a requirement that the deviation from reference temperature should be less than 12 °C. At the moment the temperatures in zone 2 and 3 are within the requirement. There are difficulties in zone 1, in section east the temperature is too high and the effects of the walking beam cause large temperature drops in section west and middle. The largest possibility for improving the temperature control will thus be in zone 1.

5.2 Controller

In this section the implementation of the different controllers whose theory was given in Chapter 3 is presented. It starts with the split range controller, then the MPC controller and finally the DMPC controller.

(48)

38 5 Control design

5.2.1 Split range controller

The split range controller has a PI-controller, which has been tuned and the version of the control signal to the burners are different than for the current con-troller. The split range controller has two different parameters, cycle time and swap time. A longer cycle time or swap time would give less switches, but at the cost of worse temperature control. Both cycle time and swap time were chosen to 45 seconds multiplied with the numbers of burners in the section which is the same as the cycle time in the current controller. For the section in which there is only one burner it is not possible to implement split range controller and in those sections the current control will remain.

5.2.2 MPC

The MPC formulation presented in Chapter 3 has been on the fundamental MPC controller with linear model, quadratic cost function and linear constraint. The flexibility of MPC allows the formulation to be modified to have a nonlinear in-ternal model, a more advanced cost function and non-linear constraints. The difference is that the optimization problem can become much more complex, but the underlying idea is still the same.

The MPC and DMPC problem were formulated by using the toolbox YALMIP [5] and was solved with the solver MOSEK [2]. To use the MPC formulation in (3.5) it needs to be some modification. The MPC controller will use the model derived in Chapter 4. The model needs to be rewritten to separate the controlled inputs, i.e. the burner from the disturbance inputs, the doors and the walking beam. The disturbance inputs are known beforehand since the walking beam is used according to a specified time and if there is a bar to be charged/discharged the doors will open before/after the walking beam is used. No constraints for temperatures in zone 1 were included in the MPC formulation since the MPC controller could not calculate control signals with the constraints.

For the current controller the control signal is converted by the burner system described in 2.2. It would be possible to let the control signal from the MPC controller be converted the same way as for the PI controllers. For the MPC controller the burners are directly controlled. The first reason is that the perfor-mance of the MPC controller would be limited by the burner system. The second reason is that the burner system would be a part of the MPC controller’s internal model of the system which would increase the complexity of the optimization problem, especially when the burner system is a time-variant system. The con-straint on the input signal needs to be modified since the burners are binary. The MPC formulation for the furnace can be written as

(49)

5.2 Controller 39 minimize u(·) N X j=1 ||_{x(k + j) − r(k + j)||}2 Q1+ ||u(k + j − 1) − u(k + j − 2)|| 2 Q2+ λ|| subject to x(k + i) = Ax(k + i) + Bu(k + i) + Bddm(k + i) xmin− ≤xi ≤xmax+  ui ∈ {0,1} 0 ≤  (5.1)

In the cost function (5.1) the term for the inputs will be zero if there were no switches during the prediction horizon. In the case of a switch for a burner the

term will always be 1 since (0 − 1)2 = (1 − 0)2 = 1. This makes it possible to set

the controller to optimize for a specific relation between a switch and deviation from the reference temperature.

Some of the parameters in (5.1) was given as xmin = r − 12, xmax = r + 12 and

λ = 109 since λ has to be large. The other parameters were decided by trying different configurations. It was done by simulating the system with the MPC controller for a half hour. This simulation was done with the same model as the internal model in the MPC controller and in the furnace. The prediction horizon

N was chosen by testing different N with Q19,9 = I

9×9 _{and Q}

217,17 = I

17×17_and

the result can be seen in

Table 5.1:The performance for the different prediction horizon.

N Mean temperature deviation Switches Max time

N = 6 57.38 179 0.25

N = 12 50.31 340 0.87

N = 18 50.85 217 7.09

N = 24 51.39 149 20.25

For N = 18 and 24 was the max time over the sample time and therefore not chosen. Between N = 6 and N = 12 was the performance better in different areas,

N = 6 had fewer switches and N = 12 better temperature control. N = 12 was

chosen since it had a better potential with a for it better ratio between Q1 and Q2.

The parameter Q2was chosen by testing different values while Q19,9 = I

9×9_{. The}

performance can be seen in Figure 5.1. It becomes a choice between switches and

mean temperature deviation and we believe that Q2 = 10 is the best choice.

It gave the following parameters N = 12, xmin= r − 12, xmax= r + 12, Q19,9 = I

9×9_,

Q217,17= 10 ∗ I

(50)

40 5 Control design

45 50 55 60 65 70 75

Mean temperature deviation from reference (celsius) 50 100 150 200 250 300 350 400 swicth Comparsion of different Q2 Q 2 = 0.1 Q₂ = 0.5 Q 2 = 1 Q₂ = 5 Q 2 = 10 Q 2 = 15

Figure 5.1:The comparison between different Q₂for MPC controller.

5.2.3 DMPC

The DMPC controller was implemented by three MPC controllers for each of the zones. The main advantage was to keep the same structure as the furnace. The burners from the same zone have a larger impact than the burners from another zone. In Figure 5.2 a sketch show how each MPC controller is working. The MPC controller for each zone will use the same MPC formulation as in (5.1). The dif-ference is for the internal model where the temperature from a nearby zone is assumed to be constant during the prediction horizon.

(51)

5.2 Controller 41

MPC 1 MPC 2 MPC 3

Zone 1 Zone 2 Zone 3

Figure 5.2:A sketch of how each MPC-controls get information.

The parameters were chosen in the same way as for the MPC controller. The result for the different prediction horizon N can be seen in the table below

Table 5.2: The performance for the different prediction horizon for DMPC

controller.

N Mean temperature deviation Switches Max time

N = 6 52.73 90 0.11

N = 12 50.65 170 0.46

N = 18 50.00 159 3.64

N = 24 49.83 168 12.54

N = 24 was not chosen since it took to long time to calculate the control

sig-nals. Between N = 6 and N = 18 the performance was better in different areas, but N = 18 was chosen since it had a better potential with more favorable ratio

between Q1and Q2.

The result for the different Q2 can been seen in the Figure 5.3. It becomes a

choice between switches and mean temperature deviation and we believe that

Q2= 30 is the best choice. The following parameters was used for the DMPC

N = 18, xmin = r − 12, xmax = r + 12, Q13,3 = I 3×3_{, zone 1:} _Q 26,6 = 30 ∗ I 6×6_, zone 2: Q25,5= 30 ∗ I 5×5_{, zone 3:} _Q 26,6= 30 ∗ I 6×6_{and λ = 10}9

(52)

42 5 Control design

50 55 60 65 70 75

Mean temperature deviation from reference (celsius) 20 40 60 80 100 120 140 160 swicth Comparsion of different Q2 Q₂ = 1 Q 2 = 5 Q 2 = 10 Q₂ = 15 Q 2 = 20 Q₂ = 25 Q₂ = 30 Q 2 = 35

Figure 5.3:The comparison between different Q2 for DMPC controller.

5.3 Simulation

The simulation was done in SIMULINK with the model based on prediction and with initial temperature, reference value and disturbance input (the doors and walking beam) from the data collection. The simulation time was one hour, which allows the temperature to weigh in. For the MPC/DMPC there were two simula-tions, one was the ideal case and the other robustness test case. In the ideal case was the internal model for the MPC/DMPC was the same as for the fur-nace. In the other case was the internal model for the MPC/DMPC estimated the same way as in Chapter 4 but the estimation data where the data from the last two thirds instead. The input temperature for the MPC/DMPC was modified by adding Gaussian noise with variance 3. The reason for it was to represent the uncertainty with the model and that the furnace is a nonlinear system. The ro-bustness test was to see if the performance changed.

The reference temperature for each section can be seen in Table 5.3 and the tem-perature from the simulation with the different controllers can be seen in Figures 5.4-5.9. More detailed picture of zone 2 and 3 for MPC and DMPC with lines where the temperature should be within are seen in Appendix B.

(53)

5.3 Simulation 43

Table 5.3:The reference temperature in Celsius for each section.

Zone Temp West Temp Mid Temp East

1 960 955 955 2 1065 1060 1060 3 1098 1098 1098 0 500 1000 1500 2000 2500 3000 3500 940 960 980 1000 1020 1040 1060 1080 1100 1120 Original Temperatures Temperature (celsius) Time (seconds) Z1 W Z1 M Z1 E Z2 W Z2 M Z2 E Z3 W Z3 M Z3 E

(54)

44 5 Control design 0 500 1000 1500 2000 2500 3000 3500 Time (seconds) 940 960 980 1000 1020 1040 1060 1080 1100 1120 Temperature (celsius) Split-range Temperatures Z1 W Z1 M Z1 E Z2 W Z2 M Z2 E Z3 W Z3 M Z3 E

Figure 5.5:The simulated temperatures with Split-range controller during 1

hour. 0 500 1000 1500 2000 2500 3000 3500 Time (seconds) 920 940 960 980 1000 1020 1040 1060 1080 1100 1120 Temperature (celsius) MPC Temperatures Z1 W Z1 M Z1 E Z2 W Z2 M Z2 E Z3 W Z3 M Z3 E

Figure 5.6: Thesimulated temperatures with MPC during 1 hour with the

(55)

5.3 Simulation 45 0 500 1000 1500 2000 2500 3000 3500 Time (seconds) 920 940 960 980 1000 1020 1040 1060 1080 1100 1120 Temperature (celsius) MPC Temperatures Z1 W Z1 M Z1 E Z2 W Z2 M Z2 E Z3 W Z3 M Z3 E

Figure 5.7: The simulated temperatures with MPC during 1 hour with the

robustness test case.

0 500 1000 1500 2000 2500 3000 3500 Time (seconds) 920 940 960 980 1000 1020 1040 1060 1080 1100 1120 Temperature (celsius) DMPC Temperatures Z1 W Z1 M Z1 E Z2 W Z2 M Z2 E Z3 W Z3 M Z3 E

Figure 5.8:The simulated temperatures with DMPC during 1 hour with the

(56)

46 5 Control design 0 500 1000 1500 2000 2500 3000 3500 Time (seconds) 920 940 960 980 1000 1020 1040 1060 1080 1100 1120 Temperature (celsius) DMPC Temperatures Z1 W Z1 M Z1 E Z2 W Z2 M Z2 E Z3 W Z3 M Z3 E

Figure 5.9:The simulated temperatures with DMPC during 1 hour with the

robustness test case.

The number of times that the burners switched between on and off during the simulation was saved. This was also done for the corresponding measurement data. It was used to compare the controllers between each other and also the controller in the actual furnace. The result can be seen in Table 5.4.

Table 5.4: The number of switches each zone had during simulation for the

different controllers and measured data for the furnace.

Controller switches z1 Switches z2 Switches z3 Total Switches

Measured 83 162 243 488 Simulated 63 133 179 375 Split range 38 67 96 201 MPC (ideal) 69 42 30 141 MPC (robustness) 172 76 64 312 DMPC (ideal) 28 11 5 44 DMPC (robustness) 66 36 64 166

Tables 5.5 - 5.7 shows the mean deviation from the reference temperature in each section for the last half hour. The last half hour was used since then the fluctuations in temperature had subsided, mostly for the simulation with current control in Figure 5.4. The tables were used to see the difference between the controllers in another way than the figures.

(57)

5.3 Simulation 47

Table 5.5:The mean deviation for the temperature in zone 1 for the last half

hour.

Controller Temp West Temp Mid Temp East

Measured 17.6 4.4 54.2 Simulated 11.3 3.4 39.8 Split range 12.0 4.1 39.7 MPC (ideal) 6.8 13.5 25.8 MPC (robustness) 6.8 10.3 29.3 DMPC (ideal) 7.0 15.8 24.9 DMPC (robustness) 8.3 8.1 32.4

hour.

Measured 2.0 1.9 6.1 Simulated 1.9 1.3 5.6 Split range 1.5 1.5 4.9 MPC (ideal) 6.4 2.8 3.2 MPC (robustness) 2.0 1.3 3.6 DMPC (ideal) 2.1 1.3 1.0 DMPC (robustness) 0.6 0.6 1

hour.

Measured 1.4 1.1 2.1 Simulated 1.1 0.9 1.2 Split range 1.2 1.2 1.5 MPC (ideal) 2.8 1.8 2.1 MPC (robustness) 1.9 0.8 3.5 DMPC (ideal) 0.3 1.6 1.5 DMPC (robustness) 3.6 2.1 0.9

Modeling and Temperature Control of an Industrial Furnace

Master of Science Thesis in Electrical Engineering

Department of Electrical Engineering, Linköping University, 2016