Injection for the Flue-gas Treatment at Fortum's Thermal Power Plant

Degree project in

Injection for the Flue-gas Treatment at Fortum's Thermal Power Plant

Joacim Sundberg

Stockholm, Sweden 2014

XR-EE-RT 2014:024 Automatic Control

Master’s Degree,

Abstract

The aim of this thesis is to investigate the possibility to improve the separation of HCl and SO 2 in the dry flue-gas treatment for boiler #3 at Fortum’s thermal power plant in H¨ ogdalen; by using a model predictive control instead of a PID controller to govern the slaked lime injection.

To achieve this an ARMAX model was derived using MATLAB’s System Iden- tification Toolbox and measurements of the incoming and outgoing levels of HCl, SO ₂ and the speed of the injection motor. The ARMAX model was then converted to a state space model which will be used as the internal model for the MPC predic- tions. The cost function was a quadratic problem which included the error between the output and the set points, the change rate of the input and the inputs deviation from a default value. The MPC uses both a feedforward and a feedback loop to es- timate the error over the prediction horizon. The controller also utilizes the ability to set constraints and tuning of the cost function weights.

In conclusion, the thesis shows that a MPC controlled lime injection is possible

and would offer some unique possibilities such as: natural constraints handling,

more intuitive live tuning for the operator and prioritized input control. However

the dry scrubber still struggles to suppress high amounts of incoming SO 2 and since

the project lacked a measuring unit for incoming HCl concentration the results

showing an improvements in pollution separation was not conclusive.

Summary

In 2011 Fortum reconstructed boiler #3 at H¨ ogdalens thermal power plant in order to produce more energy, which also lead to larger volumes of environmentally harmful flue- gases. Since the reconstruction Fortum have experienced poor separation of hydrogen chloride (HCl) and sulfur dioxide (SO 2 ) in their dry flue-gas treatment. In the dry flue-gas treatment (dry scrubber) slaked is injected into the flue-gas channel where it reacts with the HCl andSO ₂ , creating a solid substance, which can be filtered out with a filter.

This project investigated the possibility to improve the separation of HCl and SO 2

in the dry scrubber by using a model predictive controller (MPC) instead of the existing PID controller to inject lime into the flue-gas channel. The first step was to create a mathematical model which described the relationship between incoming HCl, SO 2

and slaked lime and outgoing HCl and SO ₂ . This was done by measuring the inputs and outputs and using MATLAB’s system identification toolbox to derive an ARMAX (Autoregressive Moving Average Extra signal) model. This model was then converted to a state space model to suit the purpose as an internal model for the MPC controller.

The second step was to create the MPC controller which governs the speed of the motors feeding lime into the flue-gas channel in MATLAB Simulink. The MPC controller tries to find the optimal change of input for the system in order to keep the outputs at the desired reference levels. This is done by setting up a cost function which associates a fictive cost to the deviation from the reference levels, the change of input and the deviation from a default input value. The cost function is designed as a quadratic programming (QP) problem which the MPC controller solves in order to find the optimal input to the system.

In other words the MPC controller tries to find the smallest input change which would give the smallest deviation from the desired outputs and default input. The derived controller uses both a feedforward connection and a feedback loop to estimate the sum of the anticipated and current error between the outputs and desired reference levels.

The controller also utilizes the possibility to implement constraints by limiting the range of possible inputs and change rate of the input. This means that the controller will only suggest an optimal input bound by the system constraints. The report will also present a simulation of the existing PID controller for comparison purposes.

The results presented in this report contains the derived state space model, the layout of the implemented MPC controller and the optimization code used by the con- troller along with the simulation results. It was found that the MPC controller manages marginally better control of the SO ₂ output and a more precise control of the HCl output. However these results are based on simulations and could change in a real implantation. Even with small improvements to the emission separation the MPC con- troller offers some interesting possibilities for the overall system. The MPC controller is capable of handling constraints much more natural than the PID controller. It could also offer live tuning of the controller to the operator along with prioritized input control.

The latter meaning it is possible to set different costs to multiple controllable inputs in

the cost function, which would mean that the controller can prioritize which input is the

most cost effective to change under the current circumstances.

Sammanfattning

Under 2011 utf¨ orde Fortum en ombyggnation av Panna 3 f¨ or att ¨ oka energiproduktio- nen vid H¨ ogdalens v¨ armeverk. Detta ledde till ett ¨ okat r¨ okgasfl¨ ode genom pannans r¨ okgasrening och Fortum har sedan ombyggnationen noterat en f¨ ors¨ amrad separation av v¨ ateklorid (HCl) och svaveldioxid (SO ₂ ) i den torra r¨ okgasreningen. I den torra r¨ okgas- reningen (torr skrubber) tills¨ atts sl¨ ackt kalk till r¨ okgaskanalen som sedan reagerar med HCl och SO ₂ . Reaktionen mellan kalk och HCl/SO ₂ skapar ett fast ¨ amne som kan filtreras med ett filter.

Detta projekt har unders¨ okt m¨ ojligheten att f¨ orb¨ attra separationen av HCl och SO 2 i den torr skrubbern genom att anv¨ anda en prediktive regulator f¨ or att styra kalkinmatnin- gen ist¨ allet f¨ or en PID regulator. Projektet inleddes med att skapa en matematisk modell som beskriver sambandet mellan inkommande HCl, SO 2 och kalk och utg˚ aende HCl och SO ₂ . F¨ or att ˚ astadkomma detta s˚ a m¨ attes indata och utdata f¨ or att sedan anv¨ andas i MATLAB’s System Identification Toolbox som sedan skapade en ARMAX (Autore- gressive Moving Average Extra signal) modell. Denna modell konverterades sedan till en tillst˚ andsmodell f¨ or att b¨ attre passa ¨ andam˚ alet som en intern modell i MPC regulatorn.

N¨ asta steg var att s¨ atta ihop sj¨ alva regulatorn som styr hastigheten av de motorer som matar in kalk i r¨ okgaskanalen i MATLAB Simulink. Denna regulator har till uppgift att hitta den optimala f¨ or¨ andringen av motorhastigheten som g¨ or s˚ a att utsignalen h˚ aller sig p˚ a en ¨ onskad referensniv˚ a. Detta utf¨ ors genom att st¨ alla upp en s˚ a kallad kostnads- funktion som associerar en fiktiv kostnad till att avvika fr˚ an referensniv˚ an, att f¨ oresl˚ a en stor ¨ andring av motorhastigheten eller att avvika fr˚ an en ¨ onskad motorhastighet. Kost- nadsfunktionen ¨ ar formulerad som kvadratisk problem som MPC regulatorn f¨ ors¨ oker l¨ osa f¨ or att hitta den optimala insignalen till systemet. Med andra ord s˚ a f¨ ors¨ oker reg- ulatorn att hitta den minsta ¨ andringen av motorhastigheten som bidrar till den minsta avvikelsen fr˚ an ¨ onskad motorhastighet och minsta avvikelsen mellan utsignal och refer- ensniv˚ a. Den framtagna regulatorn anv¨ ander sig av b˚ ade framkoppling och ˚ aterkoppling f¨ or att estimera summan av de nuvarande och f¨ orv¨ antade avvikelsen mellan utsignal och referensniv˚ a. Regulatorn anv¨ ander sig ocks˚ a av restriktioner som begr¨ ansar hastigheten p˚ a motorn och hur snabbt regulatorn kan ¨ andra den tidigare motorhastigheten. Detta betyder att regulatorn kommer endast att f¨ oresl˚ a en f¨ or¨ andring av hastigheten som lig- ger inom systemets restriktioner. Denna rapport kommer i j¨ amf¨ orelsesyfte ocks˚ a att presentera en simulation av den existerande PID regulatorn.

Resultaten fr˚ an denna rapport kommer att inneh˚ alla den framtagna tillst˚ andsmodellen,

en skiss ¨ over den implementerade MPC regulatorn, den kod som utf¨ or sj¨ alva optimerin-

gen samt diagram fr˚ an simuleringar av MPC och PID regulatorerna. I dessa resultat

visade det sig att MPC regulatorn lyckas ˚ astadkomma marginellt b¨ attre kontroll ¨ over

utg˚ aende SO ₂ samt en mer exakt kontroll av utg˚ aende HCl. Det skall dock noteras

att dessa resultat ¨ ar baserade p˚ a simuleringar och kan komma att ¨ andras i en verklig

implementation. ¨ Aven med sm˚ a f¨ orb¨ attringar av utsl¨ appsv¨ ardena s˚ a erbjuder MPC reg-

ulatorn n˚ agra intressanta m¨ ojligheter. En MPC regulator kan hantera restriktioner i

processen mycket mer naturligt ¨ an PID regulatorn. Den kan ocks˚ a justeras under drift

av operat¨ oren samt prioriterat val av kontrollsignal. Med prioriterat val av kontrollsignal

menas att det ¨ ar m¨ ojligt att f¨ orknippa olika kostnader till flera olika kontrollsignaler i

kostnadsfunktionen. Detta skulle medf¨ ora att regulatorn prioriterar en ¨ andringar av den

kontrollsignal som medf¨ or den minsta kostnaden under r˚ adande omst¨ andigheter.

Nomenclature

A State-space matrix

B State-space matrix

B _d Contains the second and third column of state-space matrix B, the first column have been set to zero

B _new Contains the first column of state-space matrix B, second and third column have been set to zero

C State-space matrix

D State-space matrix

D m Vector containing measured disturbances over the prediction horizon HCl _in Incoming Hydrogenchloride before the dry scrubber

HCl _out Outgoing Hydrogenchloride after the dry scrubber

H _p Prediction Horizon

H _u Control horizon

H _w States from which sample in the prediction horizon the controller penalizes devia- tions from the reference levels

I Unit Matrix

L Observer gain

Q Matrix used to increase the cost of the output error in the cost function q ⁻ⁱ Lag operator with a shift of i

R Matrix used to increase the cost of the change rate of the input in the cost function S Matrix used to increase the cost when the input is deviating from the default value

u ₀ in the cost function

SO _2,in Incoming Sulf urdioxide before the dry scrubber

SO _2,out Outgoing Sulf urdioxide after the dry scrubber

T _s Sampling time

u(k) Applied input to the system

u _o Default value on the controllable input

V (k) Cost function

x(k) State of the system at sample k y(k) Output from the system at sample k

z(k) Controllable output from the system at sample k ˆ

u(k) Calculated optimal input, not yet applied to the system ˆ

y(k) Estimated output from the system ˆ

z(k) Estimated controllable output

∆U (k) Vector containing the change of input Z Vector containing estimated outputs

Ψ Matrix used to calculate the contribution from the state x to the output z Γ Vector containing the reference levels over the prediction horizon

Θ Matrix used to calculate the contribution from change of input ∆U (k) to the output z

Υ Matrix used to calculate the contribution from previous input u(k − 1) to the output z

Ξ Matrix used to calculate the contribution from the measured disturbances D _m to the output z

E Output error

Glossary

ARX AutoregRessive eXtra signal

ARM AX AutoregRessive Moving Average eXtra signal

BJ Box–Jenkins

HCl Hydrogen chloride

ID − f an Induced Draft fan M P C Model Predictive Control

M IM O Multiple-Input and Multiple-Output

OE Output-Error

P ID Proportional, Integral and Derivative P RBS PseudoRandom Binary Sequence

P 0 The wet flue-gas treatment/wet scrubber

QP Quadratic Programming

RP M Revolutions Per Minute

RH Relative Humidity

SO ₂ Sulfur dioxide

Contents

1 Introduction 1

1.1 Explaining the Process from Boiler to Chimney . . . . 1

1.2 Problem Description . . . . 4

1.3 Goals of the Project . . . . 5

1.4 The General Structure of the Project . . . . 5

2 Designing the Model 7 2.1 Explaining the Dry Scrubber . . . . 7

2.1.1 The Chemical Reaction in the Dry Scrubber . . . . 7

2.2 Theory for Modelling the Plant . . . . 9

2.2.1 Autoregressive Moving Average Extra Signal Model . . . . 9

2.2.2 Using a State-Space Model . . . . 11

2.3 The Method for Creating the Model . . . . 12

2.3.1 Collecting Data for the Model . . . . 12

2.3.2 Creating a Model using MATLAB’s Identification Toolbox . . . . . 14

2.4 Validating the Model . . . . 17

2.5 The Model of the Plant . . . . 21

3 Designing the Controller 22 3.1 The Existing Controller . . . . 22

3.2 The Basic Idea of a New MPC Controller . . . . 23

3.3 The Theory Behind Model Predictive Control . . . . 24

3.3.1 Basic Formulation . . . . 25

3.3.2 Formulating the Cost Function . . . . 29

3.3.3 Adding Constraints to the Controller . . . . 30

3.3.4 Adding Measured Disturbances to the Prediction . . . . 31

3.3.5 Estimating the State of the System with an Observer . . . . 32

3.3.6 Finding the Optimal Solution . . . . 33

3.4 Implementing the MPC in Simulink . . . . 35

3.4.1 The First Step Towards Controlling the System . . . . 35

3.4.2 Creating the Observer . . . . 35

3.4.3 Adding the Feedforward Connection . . . . 36

3.4.4 Finding and Implementing the Constraints . . . . 37

3.4.5 Tuning the Weights in the Cost Function . . . . 39

3.4.6 The Optimization . . . . 41

3.5 Copying the Existing PID Controller . . . . 41

3.6 Simulation Results . . . . 42

3.6.1 Simulation of the MPC controller . . . . 42

3.6.2 Simulation of the PID controller . . . . 44

4.2 The MPC Controller . . . . 47

4.3 The Existing PID Controllers . . . . 48

4.4 Final Thoughts . . . . 48

5 Discussion 49 5.1 The Use of the S Matrix . . . . 49

5.2 The Installation of an HCl Measuring Unit . . . . 49

5.3 Start-up Time for the Second Feeding Line . . . . 49

5.4 Suggested Changes to the MPC . . . . 50

5.5 The Future . . . . 50

Appendices 52

A The Numerical Values for Ψ Υ Θ and Ξ 52

B The Original PID Block Diagram 55

C Optimization Code Used in the MPC 56

D A Mathematical Model of the Lime Injection 57

1 Introduction

The modern household consumes on average 12000 kWh(apartment) ¹ each year with the two largest contributors being electricity for electronics and heating of the household.

The average household also produces 465 kg of household waste ² each year and the numbers keep increasing.

1.1 Explaining the Process from Boiler to Chimney

The basic principle for boiler #3 (P3) at Fortum’s thermal power plant in H¨ ogdalen, Stockholm, is to burn household waste in order to turn water into steam. The high pressure steam is then used to drive a turbine to produce electricity. After the turbine the steam passes through a heat exchanger where it is cooled down with district heating water, the steam is now in liquid form and is fed back to the boiler for reheating. The district heating water used to cool down the steam in the heat exchanger is pumped out in order to heat residences connected to the district heating network.

We now know the ”why” and will continue with the ”how”. The waste is delivered directly by garbage trucks to the plants ”bunker”. The bunker is used as a storage station and fuel buffer for the boilers. The waste is then delivered from the bunker to the boiler inlet by a grappling claw. The inlet is shaped as a steep canal which also works as a buffer when fully loaded. At the bottom of the canal is a large hydraulic pusher which feeds waste into the boiler for combustion.

Inside the boiler the waste lies on a grate which is constantly moving back and forth, pushing the waste trough the drying zone, combustion zone, extinguish zone and finally pushing them down into a water bath where the remaining waste and ashes are extinguished and transported out of the boiler.

The warm flue-gases from the fire rises to the top of the boiler where it is led through different sections of tubes in which water is heated up and turned into steam. In order to get the flue-gases to move in the right direction and not leak out into the surrounding building, an induced draft fan (ID-fan) is placed half way down the flue-gas channel.

The ID-fan creates a suction which guides the flue-gases from the boiler, through the Economizer and dry scrubber which will be explained further down in this section.

When the flue-gases leaves the actual boiler it enters the Economizer which acts as a huge heat exchanger, cooling the flue-gases with water before they enters the dry scrubbers which is the first treatment process of the flue-gases.

In the dry scrubber slaked lime and activated carbon is fed into the flue-gas channel.

The slaked lime reacts with the acid gases, mainly sulphur oxide, SO 2 , and hydrogen chloride, HCl, while the activated carbon absorbs any heavy metals. The flue-gases then pass through a fabric filter in which particles in the flue-gases, reacted lime and activated carbon is filtered out and transported to storage silos.

1 Energir˚ adgivaren, http://www.energiradgivaren.se/2011/09/elforbrukning-i-en-genomsnittlig-villa- respektive-lagenhet/

2 Sopor.nu, http://www.sopor.nu/Rena-fakta/Sverige-jaemfoert-med-EU

The flue-gases have now reached the ID-fan and are now lead towards the chemical treatment where further separation of pollutions will take place. The flue-gases are now driven by a booster-fan which sucks the flue-gases from the ID-fan, through the chemical process and passes them to the chimney.

In the chemical cleaning process, also known as P0, the flue-gases are washed with water as they pass through. [Lundgren and Hall, 2008] states that this allows for any water-soluble substances to be separated from the flue-gases. The first step is the acid wet scrubber where any renaming HCl, hydrobromide, ammonia, hydrogen fluoride and particles are separated, hence the name acid wet scrubbed. The normal pH-level for the first scrubber lies below 1 pH. In the second scrubber, also known as the neutral scrubber, the remaining SO 2 is separated much in the same way as the first scrubber.

The flue-gases have now reached the booster-fan which passes them to the chimney

where they leave the plant. Figure 1 shows the whole process described above and

illustrates how the flue-gases move through it.

P 0

1 2

3 4

5 6

7 8 9

10 11 13 14

15 16 12

 1 . B u n ker  2 . Cra n e  3 . F eed er Ch u te  4 . F u el Pu sh er  5 . Pri ma ry Co mb u stio n A ir Fa n  6 . Gra te (In sid e Co mb u stio n Ch amb er )  7 . Pu sh er / Di sc h arg er S h af t  8 . S u p er H ea ters (S tea m Tu b es )  9 . E co n o mi zer  10 . F ab ri c F ilter / Dry Sc ru b b er  11 . S la ked L ime & A ctivated Ca rb o n S ilo  12 . F iltret Flu e- ga ses  13 . ID -F an  14 . Ch emi ca l F lu e -g as T re atmen t  15 . B o o ster -fa n  16 . Ch imn ey  Fl u e- ga s Fl o w (D irect io n al)  As h co n vey o r

 Sl ak ed Li me an d Act iv at ed C ar b o n  Pr imar y C o mb u st io n Ai r

Raw

Filtered

Clean

Figure 1: Showing the whole process from the bunker where the waste enters the plant to the chimney

where the remaining flue-gases exits the plant. The different locations throughout the process are

marked with numbers and named in the list below the figure.

1.2 Problem Description

When the household waste is incinerated it produces flue-gases with high amounts of environmentally harmful substances. In order to meet the high demands of today’s pollution restriction, these flue-gases need to be cleaned and filtered.

For boiler #3 Fortum have installed both a dry scrubber and a wet scrubber to clean the flue-gases, both which has been mentioned briefly above. In the dry scrubber hydrated lime is injected to remove the acid gases (such as SO ₂ and HCl) as well as activated carbon which absorbs any heavy metals (such as mercury among others). The flue-gases then pass through a fabric filter which filters out particles, used hydrated lime and activated carbon. The flue-gases can then continue through the wet scrubber while the waste, excess lime and carbon from the dry scrubber are taken care of.

The existing lime injection at H¨ ogdalen is done by two separate dosing lines which feed lime directly into the flue-gas channel in which lime reacts with the acid flue-gases.

The two lines consist of two sets of screws which feeds lime into transport tubes. The tubes are then pressurised by compressors which then sends the lime along the tubes into the flue-gas channel.

The way the amount of lime is regulated is by controlling the speed of the screws feeding lime into the transport tubes. This is done by a feedback loop which measures the levels of HCl and SO ₂ after the fabric filter, calculates the difference to the desired value of HCl and SO ₂ , feeds it as input to a PID regulator which controls the speed of the screws. In order to get a faster response from the regulator the HCl and SO 2 are also measured before the lime injection, and used to calculate the theoretical amount of lime needed for the reaction. This is then used to set the base speed of the screws in order to suppress the incoming amount of HCl and SO 2 . During normal levels of HCl and SO 2 , values around the desired reference level, only one of the two lines is feeding lime into the flue-gases. As soon as the levels of HCl or SO ₂ exceeds a certain threshold, the second line is fired up and the amount of lime needed is split between the two lines.

When the level of pollution drops below the threshold the second line stops and lime is once again only fed from one line.

In 2011 Fortum had Boiler #3 rebuilt and modified in order to produce more energy.

This meant an increase in the combustion which in turn led to higher volumes of flue- gases. However, the dry scrubber was rebuilt with too low capacity due to a poor estimation and is struggling to handle the large quantities of flue-gas. This has led to that the fabric filter is operating at a higher differential pressure, with poor separation of the emissions, which has resulted in a large consumption of hydrated lime at a high economic cost. The high differential pressure also affects the boilers negative pressure control which in turn gives a negative effect on the combustion and the flow of the flue- gases. A second problem is that the SO ₂ tends to peak which the current regulator handles by injecting large quantities of lime, making the regulation of SO ₂ and HCl to oscillate.

Fortum is constantly monitoring the lime consumption, the levels of HCl and SO ₂ ,

the differential pressure over the fabric filter and the overall state of the process. Ac-

cording to recent measurements the amount of lime injected into the flue-gases is double

the originally estimated levels. It is Fortum’s desire to investigate the possibility of a new controller which would reduce lime consumption while maintaining or lowering the levels of HCl and SO ₂ after the fabric filter. This would reduce the economic cost of lime and might also contribute to lowering the differential pressure over the fabric filter since less lime is being filtered.

1.3 Goals of the Project

The goal of this project will be to propose how to improve the flue-gas cleaning in the dry scrubber of Boiler #3 at Fortum’s thermal power station in H¨ ogdalen. This will be done by trying to implement a Model Predictive Control (MPC) for the lime injection.

The results expected from this project are presented below with a short explanation:

• A mathematical model which describes the relationship between incoming pollution (HCl and SO ₂ ), the amount of injected lime into the flue-gas channel and the output in the form of HCl and SO 2 after the fabric filter.

• A suggestion of a MPC which uses the above mentioned model in order to control the amount of lime injected in order to suppress levels of HCl and SO 2 below a set reference level after the fabric filter.

• A simulation of how the implementation of above mentioned MPC would function.

• In order to evaluate the efficiency of a new controller, this project will also compare the current control method with the new MPC method.

To clarify the goals of this project. The project will not cover a real implementation of a MPC method on the lime injection. This means that the MPC method will not be used on-line in the actual Dry Scrubber. All benchmark data and comparisons will be based on simulations and will suffer from any inconsistencies between the real Dry Scrubber and the simulations.

1.4 The General Structure of the Project

The project was initiated with a search for literature related to the chemistry behind the

dry scrubber and for any previous attempts to implement a MPC controller for the lime

injection. After the literature search it was decided to handle the system as a black box

and use the inputs (lime _in , HCl _in , SO _2,in ) and outputs (HCl _out , SO _2,out ) to create a

mathematical model. The reason to treat the system as a black box was because many

of the substances produced in the chemical reaction was not measured, which meant that

it would be impossible to determine how much lime that was being used in the different

reactions. Based on these initial decisions a series of experiments were conducted in order

to capture the dynamics of the system. The measured data from these experiments was

then used to create a model which could simulate and predict the output from measured

input data. The next part of the project was to design a predictive regulator which uses

the predictions from the derived model to determine the optimal control action in order

to keep the output at a set reference point. Next the MPC controller was implemented

in a simulation environment in order to evaluate its performance while controlling the

system. The final step to achieve the goals of this project was to simulate the existing

PID controller while using the same model as the MPC controller in order to compare

the two controllers.

2 Designing the Model

The following section will explain how the plant was modelled. Section 2.1 will described the functionality of the dry scrubber. Section 2.2 will give the underlying theory used in the modelling process. Section 2.3 will then explain the practical aspects of the modelling while section 2.4 will showing how the model was validated. Finally section 2.5 will present the derived state-space model.

2.1 Explaining the Dry Scrubber

When the flue-gases have left the boiler and passed through the Economizer it arrives at the dry flue-gas treatment also known as the dry scrubber. The aim for the dry scrubber is to remove the acid components and heavy metals in the flue-gases. This is done by adding slaked lime which reacts with sulphur oxide, SO ₂ , hydrogen chloride, HCl, and activated carbon which absorbs heavy metals.

The slaked lime is stored in a silo which splits into two chutes leading the lime down to a set of feeding screws, four for each line. The screws are governed by a motor and by controlling the speed of the motor one can regulate the amount of lime being fed into the flue-gases. When the screws are operating they are feeding lime into a chute with a rotary valve which feeds lime into a transport tube. The tube is pressurized by an air compressor in order to send the lime along the tube and into the flue-gas channel where it starts to react with the flue-gases.

The lime is now mixed with the flue-gases and travels with them to the fabric filter where reacted lime, activated carbon and other particles are filtered out and sticks to the surface of the filter. This creates a layer of reacted and semi-reacted lime which the flue-gases are forced to pass through due to the suction from the ID-fan. This layer actually help the flue-gases to react with the lime so that even more SO ₂ and HCl are filtered out.

In order to prevent the layer of lime and particles on the fabric filter to become too thick and hinder the flow of the flue-gases an automatic cleaning system has been installed. When the measured differential pressure over the filter rises above a certain threshold a cleaning sequence is started which releases pressurized air in short bursts down along the fabric tubes, forcing some of the layer to break off and fall down into ash conveyors.

The fabric filter is illustrated in Figure 2 and one of the two dosing lines is illustrated in Figure 3.

2.1.1 The Chemical Reaction in the Dry Scrubber

The reaction between the SO ₂ and HCl in the flue-gases and the slaked lime, Ca(OH) ₂ , is according to [Wilkstr¨ om-Blomqvist et al., 2006, GmbH, 2011] mainly described by the following two reactions:

Ca(OH) 2 + 2HCl(g) → CaCl 2 + 2H 2 O(g) (1)

Figure 2: Illustrates the dry scrubber. Slaked lime and activated carbon enters and mixes with the flue-gases right before the dry scrubber. The lime reacts with the flue- gases and sticks to the filter surface while the flue-gases passes through the tubes and up into the outlet.

Ca(OH) ₂ + SO ₂ → CaSO ₂ ∗ 1

2 H ₂ O + 1

2 H ₂ O (2)

Where the produced CaSO ₂ and CaCl ₂ are solid compounds which will be filtered out in the fabric filter.

According to [Wilkstr¨ om-Blomqvist et al., 2006] when there is an excess of calcium, reaction 1 will be replaced by:

Ca(OH) 2 + HCl(g) → CaOHCl + H 2 O (3)

The produced CaOHCl from reaction 3 will in turn start to react with the SO 2 in accordance with:

CaOHCl + SO ₂ + 1

2 H ₂ O → CaSO ₃ ∗ 1

2 H ₂ O + HCl (4)

Reaction 4 produces more HCl which will react in reaction 3, producing more CaOHCl for reaction 4 until the amount of calcium has reached normal levels at which reaction 1 once again occurs. According to [Wilkstr¨ om-Blomqvist et al., 2006] this shows that the presence of HCl can increase the level of SO ₂ separation in the dry scrubber. It is also stated in [Wilkstr¨ om-Blomqvist et al., 2006] that the produced CaCl 2 from reaction 1 increases the level of water on the surface of the lime particles which gives the SO 2

longer time to react with the lime. The CaCl ₂ also increases the pore volume of the

slaked lime which increases the contact area between the SO ₂ and slaked lime. However

Figure 3: The slaked lime is stored in a silo above the dosing screws. The amount of lime feed to the flue-gas channel is regulated by the speed of the screws. The lime passes through the rotary sluice into a pressurised tube which feeds the lime into the flue-gas channel.

it is stated in [Wilkstr¨ om-Blomqvist et al., 2006] that reaction 1 is a first order reaction and reaction 2 a second order reaction. This means that the reaction between HCl and slaked lime will take precedence over reaction 2

Another important factor for reactions 1 and 2 is the relative humidity (RH) in the flue-gases. [Wilkstr¨ om-Blomqvist et al., 2006] states that reaction 1 and 2 takes place in the complex bound water on the lime particles which means that the reaction will become more effective with a higher RH-level. However, too high RH factor will create complications in the form of corrosion and sticking tendencies in the dry scrubber. This has led to that the dry scrubber is operating at a temperature which would give the best RH-factor but still prevents moisture in the flue-gases. At Fortum H¨ ogdalen the dry scrubber is operating at 140 ^◦ C.

2.2 Theory for Modelling the Plant

This section will explain the theory behind the ARMAX and state-space model which are later used in the mathematical representation of the dry scrubber.

2.2.1 Autoregressive Moving Average Extra Signal Model

A regressive model describes the relationship between a dependent variable and a number

of independent regression variables, meaning that one can estimate a variable y from

measurements of x. For dynamical systems, y(t) is estimated from previous values of

y and u at the time t − 1, t − 2, ..., t − n. One commonly used regressive model is the

Autoregressive Moving Average Extra signal Model (ARMAX), where autoregressive means that y(t) depends on previous values of itself, y(t − 1), ..., y(t − m). Moving average means that y(k) depends on a weighted average from old values of e(t). The extra signal accounts for any measurable inputs, and creates a weighted average from u(t − 1), ..., u(t − p).

This kind of model is usually used when the system dynamic is unknown or too complex to describe physically. In that case the system will be treated as a black box and the ARMAX-model will try to describe the dynamics of the system by looking at the inputs and outputs.

A system can be described as:

y(t) = η(t) + ω(t) (5)

where η(t) is a disturbance free input to the system, ω(t) is a disturbance acting on the system and y(t) is the resulting output. If it is known that there are some dynamics acting on both input and disturbance, [Ljung and Glad, 2009] states that the system could be described as:

y(t) = G(q, θ)u(t) + H(q, θ)e(t) (6)

where

G(q, θ) = B(q)

F (q) = b 1 q ^−nk + b 2 q ^−nk−1 + · · · + b _nb q ^−nk−nb+1

1 + f ₁ q ⁻¹ + · · · + f _nf q ^−nf (7) H(q, θ) C(q)

D(q) = 1 + c ₁ q ⁻¹ + · · · + c _nc q ^−nc

1 + d ₁ q ⁻¹ + · · · + d _nd q ^−nd (8) Here q is the lag operator defined as q ⁻ⁱ u(t) = u(t−i), u(t) is the input, e(t) is considered white noise and b _nk , f _nf , c nc and d _nd are coefficients to the transfer functions.

The ARMAX model is based on (6) but with the simplification F (q) = D(q) = A(q) = 1 + a 1 q ⁻¹ + · · · + a na q ^−na , which results in the following system:

A(q)y(t) = B(q)u(t) + C(q)e(t) (9)

The system in (9) is illustrated in Figure 4 which shows that while u(t) and e(t) have their own dynamics, B(q) and C(q), they both pass through the dynamic _A(q) ¹ .

What remains now is to set the coefficients a i , b j , c _k and the structure parameters nk, na, nb, nc. The structure parameters govern how many past data samples the ARMAX model uses while the coefficients have a scaling effect on these samples. The tuning of these parameters can be done manually or in system identification programs. How to find these with the help of such programs will be explained in section 2.3.2.

(9) can now be rearranged so that it is possible to predict the output y(t) of the system.

C(q) ⁻¹ A(q)y(t) = C(q) ⁻¹ B(q)u(t) + e(t) (10)

C

∑

B 1/A

e

u y

Figure 4: Block diagram showing the ARMAX-model. The input u and noise e passes through the filters B and C separately, before passing through the filter _A ¹

y(t) =



1 − C(q) ⁻¹ A(q)



y(t) + C(q) ⁻¹ B(q)u(t) + e(t) (11) [Ljung and Glad, 2009] states that due to the properties of white noise, the term e(t) is independent of past events and cannot be determined at time t. So in order to calculate the prediction of y(t) the term e(t) is ignored in (11).

ˆ y(t) = 

1 − C(q) ⁻¹ A(q) 

y(t) + C(q) ⁻¹ B(q)u(t) (12) The right side of (12) consists of past values of y(t) and u(t), in other words the ARMAX model only uses past values of the output and input in order to predict the next output ˆ

y(t) at time t.

2.2.2 Using a State-Space Model

A common way to represent a system mathematically is to write it on state space form, see (13). The state space model will connect the input, output and state of the system through a series of differential equations.

˙x(t) = A(t)x(t) + B(t)u(t)

y(t) = C(t)x(t) + D(t)u(t) (13)

Where u(t) is the input to the system, y(t) the output and x(t) is the state of the

system. The theory, stated in [Ljung and Glad, 2009], is that if the input u(t) and the

state x(t) is known at time t = t ₀ then the change of the system ˙x(t ₀ ) is also known.

Since both the state and the change of the state is known at t = t 0 , the state of the system at t = t 0 + τ can be calculated. In [Ljung and Glad, 2009] it is also stated that the effect of all previous inputs which has placed the system in its current state are stored in the vector x(t), hence the only information needed to calculate the output is x(t 0 ) and the inputs for t > t 0 . This leads to an easy and compact way to represent the system mathematically.

2.3 The Method for Creating the Model

This section will explain the practical approach of how the system was modelled, how experiments were conducted and how the measurements were used to design the plant model.

2.3.1 Collecting Data for the Model

Fortum is constantly monitoring the plant in order to detect any abnormalities or mal- functioning equipment. They also store a huge amount of these measurements in a backlog which can be used to study trends or changes over long periods of time. This made it possible to conduct the required experiments and collect the necessary data without using any extra measuring instruments.

The main interest for this project was the incoming and outgoing amount of HCl and SO 2 along with the speed of the motors feeding lime into the flue-gas channel. However the instrument for measuring the incoming amount of HCl had been removed a while back and replaced by a constant value inside the calculations for the existing controller.

How this was handled will be discussed in Section 2.3.2.

Since the incoming amount of HCl and SO 2 could not be affected the only way of

controlling the levels of HCl and SO ₂ after the dry scrubber was to regulate the amount

of lime being fed into the flue-gases. The first step in determining how the lime affected

the emission of HCl and SO 2 was to conduct a series of step responses to determine the

time constant of the system. Only one of the two feeding lines was used during these

experiments due to the fact that both lines are identical in functionality. The assumption

here is that if one line is operating at 80% motor speed it delivers the same amount of

lime as if both lines were operating at 40% each. Another important note is that both

lines are limited to a minimum speed of 10% in order to always keep a protecting layer of

lime on the fabric filter. This means that with only one active line the minimum speed

would be 10% instead of 20% if both lines were active. The result from performing a

step from 10% to 100% in the motor speed is presented in the Figure 5.

Figure 5: Showing a step response from the system when making a step from 10% to 100% in the lime motor.

The electrical RPM given by the frequency converter governing the motor was chosen to represent the speed of the motor. The relationship between electrical RPM and percentage of motor speed is linear with 136 RPM representing 10% and 1360 RPM representing 100%. It was found that the time constant of the system lies around 5 minutes. This information was then used to conduct an experiment where the speed of the motor varied as a pseudorandom binary sequence (PRBS). According to [Ljung and Glad, 2009], if the input is changed between two levels with a random interval the input will contain a rich frequency content. It is also important to let some of the step responses die out in order to capture the system dynamics.

Two experiments were conducted in which the first used a binary signal between minimum and maximum RPM with random change interval, seen in Figure 6. The second experiment used a random change interval but with varying step heights, seen in Figure 7. Both experiments used a sampling time of 60 seconds, T _s = 60, which was set by Fortum’s analytical tool Quadra.

Figure 6: Shows the input (motor speed) for dataset 1. The step height is changed

between 10% and 100%.

Figure 7: Shows the input (motor speed) for dataset 2. The step height is varying between different levels for each change.

The data from both experiments was then imported to MATLAB’s System Identifi- cation Toolbox to be used for the identification of the model, see Section 2.3.2.

2.3.2 Creating a Model using MATLAB’s Identification Toolbox

It was mentioned in Section 2.3.1 that the instrument measuring the amount of HCl before the filter, the incoming HCl, had been previously removed. In order for the identification toolbox to recognise that there should be an input containing information about the HCl concentration before the dry scrubber, the data had to be reconstructed.

By comparing the amounts of HCl in the raw flue-gases from the other boilers it was estimated that the levels of HCl in the raw flue-gases for boiler #3 should be around 450mg/N m ³ . To represent this a random signal was created with a mean value of 450 and a standard deviation of 10. This is a rough estimation of the actual HCl concentra- tion after the boiler and does not represent any fluctuations which would otherwise be expected. The risk with trying to add such fluctuations would be that the identification toolbox might find incorrect correlations which could affect the overall model. It should be noted that the manufactured HCl input will have a negligible impact on the final model but will visualize how a real situation might look like.

The data used in the system identification toolbox is on the form 3 inputs (RP M motor ,

HCl _in , SO _2,in ) and 2 outputs (HCl _out , SO _2,out ) with a sampling time T _s = 60 sec. The

two data sets can be found in Figures 8 and 9.

(a) Showing dataset 1 from HCl perspective

(b) Showing dataset 1 from SO 2 perspective

Figure 8: Shows the inputs and outputs for dataset 1. Figure 8a shows the HCl input and output along with the motor speed. Figure 8b shows the SO 2 input and output along with the motor speed

(a) Showing dataset 2 from HCl perspective

(b) Showing dataset 2 from SO ₂ perspective

Figure 9: Shows the inputs and outputs for dataset 2. Figure 9a shows the HCl input

and output along with the motor speed. Figure 9b shows the SO ₂ input and output

along with the motor speed

Once the data had been inserted into the System Identification Toolbox the mean values of both sets were removed, centering the data around zero. The toolbox also presents the option to merge data in order to use multiple measurements for one model.

By merging dataset 1 and 2 the toolbox actually uses both when fitting a model to the system. The next step in the modelling process was to decide the form and complexity of the model. A rule of thumbs is to always start with the simplest model and move on up until a desired accuracy has been reached. A third dataset (dataset 3) was measured and imported for validation purposes, see Figure 10, which uses a random change interval along with a varying step height.

(a) Showing dataset 3 from HCl perspective

(b) Showing dataset 3 from SO 2 perspective

Figure 10: Shows the inputs and outputs for dataset 3 which is the verification set.

Figure 10a shows the HCl input and output along with the motor speed. Figure 10b shows the SO 2 input and output along with the motor speed

After testing a few models, such as ARX, OE and BJ, it was found that the ARMAX- model (described in Section 2.2.1) represented the system the best. Not only did it manage to get a good fit on dataset 3 after the mean values had been removed, but it also maintained a good fit even without removing the mean values. Something the other models failed to do. The validation of the ARMAX-model on dataset 3 can be found in section 2.4.

MATLAB can now be used to convert the ARMAX model into an observable state

space model, see Section 2.2.2 for explanation. The state space form is much easier to

handle then the ARMAX representation and by making it observable it is possible to

use an observer to estimate the state of the model later on. The state space model can

be found under the results in Section 2.5.

2.4 Validating the Model

Once the model had been completed it was time to evaluate how well the model represents the real system. The plots in Figure 11 and 12 shows the measured and simulated output with their respective ”fit” value. This value acts as a measurement of how well the simulated output corresponds to the real output of the used dataset.

It can be found that the HCl model has a poor fit to the real HCl output. Since the HCl input only consists of an estimated mean value and some added noise, it is impossible for the system identification toolbox to find any correlation between the HCl input and output. Instead the toolbox uses more of the SO 2 input and the speed of the lime motor to estimate the HCl output. However it can be seen in Figure 13 that the model simulates an increase in HCl _out when the motor speed drops and a decrease of HCl out when the motor speed increases. This shows that the HCl model has at least captured the basic trends of the outgoing HCl.

The SO ₂ model on the other hand has a quite good fit in both cases. This is due to the fact that the identification toolbox is able to use information from both the incoming and outgoing SO ₂ as well as the speed of the lime motor in the modelling process.

(a) Showing the simulated and measured HCl output with mean values removed along with associated fit value.

(b) Showing the simulated and measured SO 2 output with mean values removed along with the associated fit value.

Figure 11: Showing the simulated and measured output of both HCl and SO ₂ while

using dataset 3 without mean values. The figure also shows the corresponding fit values

which acts as a measurement of how well the model represented the real system.

(a) Showing the simulated and measured HCl output without modifications to the dataset along with the associated fit value.

(b) Showing the simulated and measured SO 2 output without modifications to the dataset along with the associated fit value.

Figure 12: Showing the simulated and measured output of both HCl and SO ₂ while using dataset 3 with mean values still present. The figure also shows the corresponding fit values which acts as a measurement of how well the model represented the real system.

Figure 13: Shows the simulated and measured HCl output along with the speed of the lime motor while using dataset 3. It can be seen that the model simulates an increase in the simulated HCl output when the motor speed drops and vice versa.

Another way of validating the model is to look at its residuals. These are basically the difference between the measured and predicted output, also known as the prediction error. For a good model the residuals should be small and uncorrelated with the inputs.

If however the residuals are correlated with the input it could mean that parts of ˆ y(k) has

not been properly modelled. Similarly if there is any correlation between the residuals

themselves it indicates that ˆ y(k) could have been better predicted. The residuals for the

derived model can be viewed in Figure 14 along with a 99 % confidence interval.

(a) Showing the residuals motor speed to HCl out. (b) Showing the residuals motor speed to SO 2 out.

(c) Showing the residuals HCl in to HCl out. (d) Showing the residuals HCl in to SO 2 out.

(e) Showing the residuals SO 2 in to HCl out. (f) Showing the residuals SO 2 in to SO 2 out.

Figure 14: Showing the residuals for the model.

2.5 The Model of the Plant

When the ARMAX model had been validated it was converted to a state-space model. The final model used in the MPC controller and to simulate the plant can be found below:

A =

























0.1880 0.5754 0 0 0 0 0 0

−0.5754 0.1880 0 0 0 0 0 0

0 0 0.3250 0.5678 0 0 0 0

0 0 −0.5678 0.3250 0 0 0 0

0 0 0 0 0.9194 0 0 0

0 0 0 0 0 0.5695 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

























B =

























0.0089 0.0147 −0.1633 0.0001 0.0281 −0.1637 0.0156 0.0376 −0.3555

−0.0387 −0.0535 0.6500

−0.0089 0.0122 0.0071

−0.0194 −0.0457 0.4360

0 0 0

0 0 0

























C =  0.5261 −0.0991 0.172 0.1866 0.7801 0.0783 −0.5145 0.0341

−0.4056 −0.4410 0.6113 0.1541 0.2323 0.7929 0.2543 −0.5902



D = 0 0 0 0 0 0



(14)

Substance Molar mass [g/mol]

Ca(OH) ₂ 74.093

HCl 36.46

SO 2 64.066

Table 1: The molar mass for Ca(OH) ₂ , HCl and SO ₂ .

3 Designing the Controller

The following sections will explain how the controller was design. Section 3.1 will describe how the existing controller operates. Section 3.2 will explain the main idea with the new controller and how it operates. Section 3.3 will summarize the theory behind model predictive control and section 3.4 will described how the MPC was implemented in Simulink.

3.1 The Existing Controller

The current controller for injecting lime into the flue-gas channel consists of two PID controllers (one for HCl and one for SO ₂ ) which uses both feedforward and a feedback loop to regulate the level of HCl and SO 2 after the dry scrubber. The feedforward is done by measuring the level of HCl and SO ₂ before the dry scrubber and adjusting them for the flue-gas flow. The theoretical amount of lime needed for reactions 1 and 2 is now calculated from (15) and (16) and added together.

Required − lime _HCl = HCl in

2 ∗ M _HCl ∗ M _{Ca(OH) 2} (15) Required − lime SO 2 = SO 2,in

M _{SO 2} ∗ M _{Ca(OH) 2} (16)

Where M _{Ca(OH) 2} , M _HCl and M SO 2 is the corresponding molar mass of the substance found in Table 1. Note that the multiplication with 2 in the denominator of (15) comes from the fact that two moles of HCl reacts with one mole of Ca(OH) 2 , which is in accordance with reaction 1.

The feedback loop measures both HCl and SO ₂ after the fabric filter and adjusts

them to 11% O ₂ and the current level of H ₂ O. The O ₂ adjustment is a standard set by

the government. Many of the emissions are measured by comparing them to a reference

gas which often contains 11% O ₂ . This also enables comparisons between different plants

and the use of standardized charts. The same can be said for the H ₂ O adjustment and

the use of N m ³ (normalized cubic meter) instead regular m ³ . The normalized cubic

meter has been adjusted for the current pressure and temperature. The adjusted values

are then fed to two PID controllers, one for HCl and one for SO ₂ . A switch then

compares and chooses which controller requires the largest change (represented by a

scaling factor between 0.5 - 4). The chosen factor is then multiplied with the required

lime from the feedforward calculations and sent to the motors driving the feeding screws.

In order to set the motor speed a simple mapping between required lime and motor speed is performed. Only one of the two feeding lines will be operating as long as the amount of required lime is less than 300 kg/h but as soon as more lime is needed the second line will start up and the required amount is split between the two lines. Both lines will be operating for a minimum of 15 minutes before going back to one line or if the required amount of lime goes below 150 kg/h.

Another important note is that the parameters of the PID controllers will change if the wet flue-gas treatment, P0, is being by-passed, for example during renovations. If this is the case the controllers will become more aggressive and the reference levels for HCl and SO 2 will be lowered in order to compensate for the lack of a second flue-gas treatment. The existing controller is illustrated in a block diagram in Figure 15.

Figure 15: Illustrates the existing PID controller. The controller uses both feedforward and a feedback loop where the feedforward makes a rough calculation of required lime and the feedback loop makes corrections depending on the output. The setpoints and controller parameters are changed if the wet flue-gas treatment is bypassed.

3.2 The Basic Idea of a New MPC Controller

Before presenting the theory and the actual implementation of the MPC this section will explain the basic idea of the new controller and how it operates.

The basic control problem is to control the levels of outgoing HCl and SO ₂ by

changing how much lime that is being added to the flue-gases. The new controller will

use the derived model from section 2, along with measurements of the current speed

of the lime motor and the incoming amount of HCl and SO ₂ , in order to predict the

resulting amount of outgoing HCl and SO 2 . For this purpose the incoming amount of

HCl and SO 2 , and the speed of the lime motor will be viewed as the input to the system.

The outgoing amounts of HCl and SO 2 will be viewed as the output. Since the speed of the lime motor is the only controllable input, the controller will now use the output prediction to find the best control action, i.e. how to change the speed of the lime motor to keep the output below the set reference levels. In this project the aim is to keep the outgoing HCl level below 80 mg/N m ³ and the SO 2 below 250 mg/N m ³ . In order for the controller to determine the most optimal control action it first needs something to evaluate the chosen action with. This is done by creating something called the cost function which basically is a summation of the output error and the change rate of the controllable input. In this project the cost function will consist of how much the HCl and SO 2 deviates from their reference values and how much the speed of the lime motor has changed. This means that when the MPC suggests a change to the input it uses the plant model to calculate how this would affect the output, and by inserting these calculations into the cost function it can also measure the efficiency of the suggested action. In other words the cost function will increase if the controller suggests a large change of the motor speed or if this change leads to a large deviation between the output and the reference levels. The controller will therefore select the change of input which leads to the smallest cost. In addition to the cost function the MPC controller will also take any constraints of the system into account. In this project the speed of the lime motors is bound to operate between 10% and 100%. It was also decided to limit the change rate of the motors to a maximum of 1000 RPM per change. The implication of this is that while the MPC is searching for the change of input which gives the smallest cost function it will never suggest an input which violates these constraints.

This can be summarized as the controller is searching for the smallest input change (change of the lime motor speed) which minimizes the output (HCl and SO ₂ ) deviations from their corresponding reference values over the prediction horizon, while at the same time maintaining certain constrains on the system.

A more mathematical description would be:

∆M otorSpeed min kr(k + i) − ˆ z(k + i)k

| {z }

Output error

+

Change of motor speed

z }| {

k∆ˆ u(k + i)k

subject to :

136rpm ≤ M otorSpeed ≤ 2720rpm

− 1000rpm ≤ ∆M otorSpeed ≤ 1000rpm

(17)

3.3 The Theory Behind Model Predictive Control

The PID controller is the most commonly used controller due to its simple design and

implementation. The basic way a PID controller operates is to measure the output from

a system, a valves position or a level in a tank, and compare it to the set point of the

system in order to calculate the error of the system, error = setpoint − output. The

error is then fed to the controller which adjusts the input in order to bring the system

towards the set point. The output is then measured again and the process is repeated.

The PID controller only answers one question: How should the input be changed in order to steer the system towards the set point at this exact time. The controller lacks any knowledge of the system it controls which means that even if a valve has reached its maximum position the PID controller will still tell it to go further as long as the error is non-zero.

As an alternative to the PID controller, the Model Predictive Control (MPC) method has become more frequently used as the computational power has increased. The MPC controller uses a mathematical model over the plant in order to predict the effect of the control action over a finite time horizon. It then uses this information to determine the most effective control action for the system. In [Maciejowski, 2002] it is stated that not only does the MPC controller handle multiple input/multiple output (MIMO) systems naturally, it can also be programmed to take account for any constraints in the system.

The following sections will explain the theory behind the MPC controller and section 3.4 will explain more about the actual implementation of an MPC controller.

3.3.1 Basic Formulation

In order for the MPC to calculate the optimal control action it first measures the output of the system, y(k), at time sample k. It then, just as the PID controller, calculates the difference between the output and the desired reference point of the system. Once the error is known it uses a model over the system to see how different changes to the input,

∆ˆ u(k), would affect the output of the system. Here the ∆ sign represent the difference or change of ˆ u(k) and the circumflex over u(k) indicates that it is the suggested optimal value which hasn’t been applied to the system yet. Another important thing to note is that the MPC doesn’t just predict the next output, y(k + 1), it actually predicts the effects of ∆ˆ u(k) over the whole prediction horizon ˆ y(k + 1), . . . , ˆ y(k + H p ). According to [Maciejowski, 2002], when the optimal change of the input has been found it is applied to the previous input, seen in (18), and u(k) is then used as input to the system.

u(k) = u(k − 1) + ∆ˆ u(k) (18)

The MPC now measures the output y(k +1), which is the result of the applied control action, and the whole process is repeated. When the controller is trying to find the next optimal solution it still takes into account the following H p samples meaning that the horizon is constantly moving ahead of the optimization. This is called ”receding horizon”.

The effect of this is that the MPC is constantly trying to predict how the system will respond to suggested inputs over the prediction horizon, and then use this information to find the optimal control action in order to keep the system at the reference trajectory.

So how does the MPC predict the output over the prediction horizon? If a system is represented by the state space model in (13), [Maciejowski, 2002] states that the future values can be found by iterating through (13). For the moment, it will be assumed that there are no disturbances and that the state can be measured.

The change of state is calculated by:

ˆ

x(k + 1) = Ax(k) + B ˆ u(k) (19)

The change of state two samples into the future would then be:

ˆ

x(k + 2) = Aˆ x(k + 1) + B ˆ u(k + 1) (20) and if ˆ x(k + 1) in (20) is replaced by (19) the resulting state change would become:

ˆ

x(k + 2) = A ² x(k) + AB ˆ ˆ u(k) + B ˆ u(k + 1) (21) [Maciejowski, 2002] also states that it is possible to use this iteration method throughout the whole prediction horizon which leads to the following equation:

ˆ

x(k + H _p ) = A ^{H p} x(k) + A ˆ ^{H p −1} B ˆ u(k) + · · · + B ˆ u(k + H _p − 1) (22) It was briefly mentioned that the MPC controller doesn’t calculate the actual input but rather how the input should change. So in order for the MPC to predict how the change of input affects the system, (19) through (21) must be expressed in terms of

∆ˆ u(k) rather then ˆ u(k). According to [Maciejowski, 2002] this can be done by first expressing ˆ u(k) in terms of ∆ˆ u(k) and previous inputs.

ˆ

u(k) = ∆ˆ u(k) + u(k − 1) ˆ

u(k + 1) = ∆ˆ u(k + 1) + ∆ˆ u(k) + u(k − 1) .. .

ˆ

u(k + H _u − 1) = ∆ˆ u(k + H _u − 1) + · · · + ∆ˆ u(k) + u(k − 1)

(23)

Here H u is introduced as the control horizon. The MPC can actually optimize different inputs to be applied at time samples k, k + 1, . . . , k + H _u . However it is common to only apply ˆ u(k) to the system and then repeat the optimization for ˆ u(k + 1).

If the calculations from (23) are applied to (19)-(21), the following predictions are

obtained:

ˆ

x(k + 1) = Ax(k) + B[∆ˆ u(k) + u(k − 1)]

ˆ

x(k + 2) = A ² x(k) + AB[∆ˆ u(k) + u(k − 1)]

+ B[∆ˆ u(k + 1) + ∆ˆ u(k) + u(k − 1)]

= A ² x(k) + (A + I)B∆ˆ u(k) + B∆ˆ u(k + 1) + (A + I)Bu(k − 1) .. .

ˆ

x(k + H _u ) = A ^{H u} x(k) + (A ^{H u −1} + · · · + A + I)B∆ˆ u(k)

· · · + B∆ˆ u(k + H _u − 1) + (A ^Hu−1 + · · · + A + I)Bu(k − 1) ˆ

x(k + H _u + 1) = A ^{H u +1} x(k) + (A ^{H u} + · · · + A + I)B∆ˆ u

· · · + (A + I)B∆ˆ u(k + H _u − 1) + (A ^{H u} + · · · + A + I)Bu(k − 1) .. .

ˆ

x(k + H _p ) = A ^{H p} x(k) + (A ^{H p −1} + · · · + A + I)B∆ˆ u(k)

· · · + (A ^{H p −H u} + · · · + A + I)B∆ˆ u(k + H _u − 1) + (A ^{H p −1} + · · · + A + I)Bu(k − 1)

(24)

The predictions in (24) can be expressed in matrix form, see (25), which is an easier format to handle.



















 ˆ x(k + 1)

.. . ˆ

x(k + H _u ) ˆ

x(k + H u + 1) .. . ˆ

x(k + H p )





















=



















 A

.. . A ^{H u} A ^{H u +1}

.. . A ^{H p}



















 x(k) +



















 B

.. . P H u −1

i=0 A ⁱ B P H u

i=0 A ⁱ B .. . P H p −1

i=0 A ⁱ B





















u(k − 1)+



























B . . . 0

AB + B . . . 0

.. . . .. .. . P H u −1

i=0 A ⁱ B . . . B P H u

i=0 A ⁱ B . . . AB + B .. . .. . .. . P H p −1

i=0 A ⁱ B . . . P H p −H u

i=0 A ⁱ B

































∆ˆ u(k) .. .

∆ˆ u(k + H u − 1)







(25)

Equation system (25) shows how to predict the state change. The predictions of the

output are now obtained, according to [Maciejowski, 2002], by iterating the second line

in (13).