Active Sensing Test Cases in Industrial Processes

ACTIVE SENSING TEST CASES IN INDUSTRIAL PROCESSES

Faculty of Engineering and Natural Sciences (ENS) Master’s thesis Examiners: Prof. Risto Ritala, Prof. Matti Vilkko September 2020

Mikko Salo: Active sensing test cases in industrial processes Master’s thesis

Tampere University

Master’s Degree Programme in Automation Engineering September 2020

The effect of optimizing the information state of an LQG-control problem was stud- ied. The measurements of the processes concerned were limited but controllable and the performance of the process was tested with different choices of measurements.

Two test cases were considered: froth flotation process and paper machine basis weight control. In the froth flotation process, a linear process model was identified and a mathematically optimal solution for the measurement optimization problem was derived and tested with the model. The effect of the choice of measurement scheduling was tested and the cost of operation was identified to be in the scope from lower limit of 0.48 with standard error of 3.63∗10⁻⁴ to upper limit of 0.74 with standard error of 3.96∗ 10⁻⁴ (N = 10⁶) with different choices of measure- ments. In the case on paper machine basis weight control, a simulator was used and a scheduling algorithm for a moving scanner measurement was built to dampen the effects of temporal and spatial deviations being intermingled due to the motion of the scanner measurement. Temporal oscillations in the basis weight detected as spatial oscillations in the estimates were greatly reduced due to the scheduling.

Keywords: LQG-control, measurement optimization, separation principle, froth flotation, basis weight control

The originality of this thesis has been checked using the Turnitin Originality Check service.

Mikko Salo: Aktiiviaistin testitapauksia teollisuuden prosesseissa Diplomityö

Tampereen yliopisto

Automaatiotekniikan maisteriohjelma Syyskuu 2020

Informaatiotilan optimoinnin vaikutusta LQG-ongelmissa tutkittiin. Tutkittavien prosessien mittausvaihtoehdot olivat rajalliset mutta valittavissa. Prosessien suori- tuskykyä testattiin erilaisilla mittausten valinnoilla. Kahta testitapausta tutkit- tiin: mineraalien rikastusprosessi ja paperikoneen neliömassan säätöprosessi. Mi- neraalien rikastusprosessissa tunnistettiin lineaarinen prosessimalli ja matemaat- tisesti optimaalinen ratkaisu johdettiin mittausoptimointiongelmalle, ja ratkaisua testattiin lineaarisella prosessimallilla. Mittausten valinnan vaikutusta testattiin ja määritellyn operoinnin kustannuksen alarajaksi saatiin0.48keskihajonnan keskiar- volla3.63∗10⁻⁴ja ylärajaksi0.74keskiarvon keskihajonnalla3.96∗10⁻⁴(N = 10⁶) eri mittausjonojen valinnoilla. Paperikoneen neliömassan tapauksessa käytettiin simu- laattoria ja kehitettiin liikkuvalle skannerille aikataulutusalgoritmi, jolla vaimen- nettiin skannerin liikkeestä johtuvaa spatiaalisten ja ajallisten häiriöiden sekoit- tumista. Neliömassassa spatiaalisena havaitut ajalliset häiriöt saatiin vaimennettua pois tilaestimaateista skeduloinnin avulla.

Avainsanat: LQG-säätö, mittausoptimointi, separaatioperiaate, mineraalien rikas- tus, neliömassan säätö

1 Introduction . . . 1

2 LQG-control . . . 3

2.1 Process model . . . 3

2.2 Kalman filter . . . 4

2.3 Cost function . . . 5

2.4 Separation of the cost function . . . 6

3 Test cases . . . 12

3.1 Froth flotation process . . . 12

3.1.1 HSC-simulator . . . 12

3.1.2 Process and control . . . 12

3.1.3 Measurements . . . 14

3.1.4 Test case objective . . . 14

3.2 Paper machine basis weight control . . . 15

3.2.1 Process and control . . . 15

3.2.2 Measurements . . . 16

3.2.3 Process model . . . 17

3.2.4 Process simulator . . . 20

3.2.5 Test case objective . . . 21

4 Application in froth flotation . . . 22

4.1 Linear process model . . . 22

4.2 Model parameters . . . 25

4.3 Tests and results . . . 26

4.3.1 Test environment . . . 26

4.3.2 Linear model accuracy . . . 27

4.3.3 Weighting factor effects . . . 29

4.3.4 The effect of measurement optimization . . . 32

5 Application in basis weight control . . . 34

5.1 Problematic disturbances . . . 34

5.2 Tests and results . . . 37

5.2.1 Dampening of MD oscillations . . . 37

5.2.2 Distinguishing between MD and CD disturbances . . . 39

6 Conclusions and future work . . . 40

References . . . 41

7 Appendix A: Froth flotation model dynamics matrices . . . 42

1 Introduction

Any real life process can be thought to be consisting of two parts: the part how the state of the process actually is, called the real-life state, and the part of how the operator thinks the state of the process is, called the information state. There is always uncertainty about the real-life state of the process due to incomplete obser- vations or uncertainty of the observations. This uncertainty is taken into account in information state. In decision making, such as feedback process control, we need to make decisions on how to control the real-life state of the process according to the information state available to us. For example when driving a car, one makes a decision to change the lane after looking into the rear and side view mirrors. Even though a car is not visible behind them, there may be a smudge in the mirror hiding one, or a car might appear after the driver ceases to observe the mirrors. For this reason the quality of the information state is equally important to obtaining high performance process control as the actual control law used. In this thesis we look into ways of optimizing the quality of the information state in industrial processes that already have been implemented with feedback control acting based on the in- formation state, usually considering the real-life state to be the most likely state according to the information state. These processes have incomplete measurement options, meaning that not all of the states of the processes can be measured at the same instance of time. Hence any unmeasured state at a given time might change due to process noise or external disturbances, and the real-life state and the most likely state according to the information state will very likely differ, causing sub- optimal control actions on the process. This leads to the question: which states should we measure at a given time instance and information state to obtain the best performance for the process?

In this thesis we will review a solution for this question in the case of an LQG- control problem. In this case, there exists an exact solution for the question. This means that the solution will give the best performance on the average for opti- mal feedback control acting based on the information state. LQG-processes are stochastic, meaning there exists random variables that alter the real-life state and information state of the process, hence the solution to the control problem will be to optimize the expected value of the quality of the control. The LQG-control problem will be formulated and solved in Chapter 2. Similar approach to the LQG measure- ment optimization has been studied by for example Meier, Peschon, and Dressler 1967; Suominen 2011. The emphasis will be on optimizing the information state, not so much on the actual feedback control of the process.

The solution to the information state optimization problem will be tested in two

use cases. In one case, a froth flotation process is analyzed as a state-space model and then formulated into an LQG-control problem. The goal is to produce as much as possible of the enriched end product, while keeping its purity to standards. The process will consist of feedback control with three possible qualities to measure, but only one device to measure them. An optimal sequence of measurements that minimizes the deviation from desired operating point is then derived, assuming the optimal control law is obeyed. The effects of the solution will be tested. The use case of froth flotation process is introduced in Section 3.1 and analyzed in Chapter 4.

The second case is the basis weight control of a paper machine. The goal is to regulate the basis weight of the paper to given standards across the wide paper web. The problem is that the scanner measuring the web basis weight can only measure a very small subset of the width of the web. The question becomes how the scanner should be moved to obtain a high quality information state across the entire web and to avoid pitfalls caused by some problematic external disturbances typical to the process. A solution to one of the worst case scenario disturbances will be derived. This disturbance is such that it makes the actuators work against the goal of the process by making a temporal disturbance appear as a spatial disturbance in the primary scanner measurement. A scheduling algorithm of the scanner will be derived that will dampen the effects of this disturbance. This use case will be introduced in Section 3.2 and analyzed in Chapter 5.

2 LQG-control

In this chapter we will formulate the Linear Quadratic Gaussian (LQG) control prob- lem and prove the optimal solution for it. Like the name suggests, LQG-problem consists of a linear process model affected by Gaussian white noise (see Insua, Rug- geri, and Wiper 2012, p. 11), and the objective of the control is to minimize a cost function quadratic in the deviation from a desired operating point. In practice, most processes are non-linear, and the LQG-control is usually applied to a linearized model of the process around some nominal point of operation, usually the desired operating point. The goal of the LQG-control will then be to regulate the process to remain as close as possible to the point of operation by applying as little control as possible. The metric defined by the cost function defines what counts as close to the point of operation and what amount of control is significant. Deviations from the point of operation can be caused by process or measurement noise, model inaccuracy, or external disturbances and the control is used to counteract these deviations. LQG-control problems can be continuous or discrete in time. Usually discretization is required due to the fact that the actual process is continuous, while the measurement and control devices are digital (Franklin, Powell, and Workman 1998, p. 57).

2.1 Process model

The system dynamic equations are given as a discrete state space model x_n+1 =Ax_n+Bu_n+ϵ_n ϵ∼N(0,Σ^(p))

y_n =C_nx_n+ν_n ν ∼N(0,Σ^(m)n ). (2.1) The vectors x_n ∈R^k, where k is the number of states, represent the real-life states at each time instant n. In LQG-problems the states usually represent deviations from a nominal point rather than absolute values. States might be actual physical states like velocity or position, or pseudo-states whose linear combinations make a physical output signal. The states are unambiguous but are not known to a certainty.

Vectorsu_n∈R^m are the control signals, wheremis the number of controls available.

Vectorsy_n ∈R^ln are the measured outputs, andl_n< k∀nis the number of measured outputs depending on the choice on measurementC_n ∈R^ln×k. The random vectors ϵandνare the process and measurement noises, respectively. In LQG-problems, the noises are zero-mean Gaussian white noise with covariance matrices Σ^(p) and Σ^(m)n

that are usually diagonal. This means that their expected values are zero and that the individual values of the noise component at each time step are independent of

its values in the past or in the future, and that there are no cross-correlations in the noises between the different states or measurements (See Åström and Wittenmark 1997, p. 380). The measurement covariance matrix depends on time, since different measurements have different covariance matrices.

The matrices A∈R^k×k andB ∈R^k×m define the system dynamics: for example how the states evolve wrt. control (B) or wrt. previous values of the states (A).

Note that there may exists multiple measurement matrices [C_n], one of which is chosen for any given given timen. By the choice of measurement at each time step (if possible), we can affect the outcome of the states and the efficiency of the control due to more accurate estimates of the states.

2.2 Kalman filter

In practice, the states cannot be known exactly due to for example process noise, measurement noise, incompleteness of measurements or unknown external distur- bances. We need an estimate of the states in order to even attempt to control them.

Such an estimate is provided by a Kalman filter. Kalman filter determines how the the states develop due to process and control dynamics, and due to measurements obtained. If we assume that states are normally distributed, the process is linear and the process and measurement noises are white Gaussian noise with zero mean, then the information about the states will continue to be normally distributed as the system evolves in time. We will now state without proof how the estimate is changed due to process dynamics and measurement, as the Kalman filter is very broadly researched (For example Burl 1999; Glad and Ljung 2000). The steps of Kalman filter are called predict and update. These are both a direct result of apply- ing Bayes’ theorem (see Insua, Ruggeri, and Wiper 2012, p. 11) and properties of a normally distributed random vector. Let us assume that the process model is of the form 2.1 and that the information about the state at current time n is distributed

x_n ∼N(µ_n+,Σ_n+).

The plus sign in the time subscript indicates that a measurement has been made at time n, and the mean and variance have been updated with measurement data obtained at that time instant. Similarly we will use a minus sign in the time subscript for mean and variance that have not yet been updated by a measurement done in that time step. The predict step of the Kalman filter is then

µn+1− =Aµn++Bun

Σ_n+1− =AΣ_n+A^⊺+ Σ^(p). (2.2)

As indicated by the minus signs in the subscripts, the predict step gives the mean and variance of the state information for the next time step without any measure- ment from the actual process at that time step. The problem with just predicting is that the process contains a noise term, which will gradually increase the uncertainty of the estimates, even if the model would be perfectly accurate and there were no other external disturbances. This can be helped by measuring the real system and updating the state information according to the measurement result. If a measure- mentC_n+1 is used, a result y_n+1 is obtained and the measurement has uncertainty dictated by the measurement covariance matrixΣ^(m)_n+1, then the state information is updated according to the following Kalman filter update step. This is again a result of Bayes’ theorem and the properties of normal distribution.

L_n+1Σ_n+1−C_n+1 C_n+1Σ_n+1−C_n+1^⊺ + Σ^(m)_n+1−1

µ_n+1+ =µ_n+1−+L_n+1 y_n+1−C_n+1µ_n+1− Σn+1+ = I−Ln+1Cn+1

Σn+1−

(2.3)

Here µ_n+1+ is the estimate mean that will be used to obtain the next time step control signal.

2.3 Cost function

Equation 2.1 defines how the states and outputs evolve in time and due to control signals and measurements. In order to have a control problem, we need to specify a goal. For a discrete-time LQG-problem, the goal is in the form of a quadratic finite horizon cost function

V_N[µ_n+,Σ_n+] = minE (_N−1

X

i=0

(x^⊺_n+i+1Qx_n+i+1+u^⊺_n+iRu_n+i) )

, (2.4)

where Q and R are positive semi-definite symmetric matrices. The elements of Q define how much deviations in the statesx_n add to the cost function. Similarly the elements of R penalize the control signals used. The optimization horizon N will define how far into the future the states and control are considered in the cost func- tion. With a largerN, the states and control will be predicted further in time, but the calculations will grow more complex also. The cost function is also a function of the current estimate and its uncertainty, since it is used to calculate the optimal control action and measurement again at each time step based on current informa- tion state. Even though the cost function includes control and measurement choices into the future up to the horizonN, only the first actions are implemented and the optimal actions are recalculated at the next time step based on new measurement data. The expected value in the cost function must be taken, since the evolution of

the states is affected by noises that are random vectors.

The matricesQandRare tuning parameters. The control weightingRis usually set small in comparison to the states weighting, but it is necessary for example so that the control signals do not grow too large and cause actuator saturation or cause the process to drift too far out of the point of linearization, where the model will no longer hold. Alternatively, in some processes the control will actually be expensive also, and it is natural to penalize the use of it. As an example of this would be a combustion engine with fuel flow as control. The tuning parameters Q and R are usually set by someone working to benefit from the process since they can tell which outputs are important and which are not. For this reason it is usually helpful to give the cost function as a function of the measured outputs rather than the states. This is because sometimes the states are not very intuitive and for the decision makers it is easier to tell what they want to emphasize. The cost function in terms of the data obtained from all the measurable outputs is given as

V_N = minE (_N−1

X

i=0

(y^⊺_n+i+1Q_yy_n+i+1+u^⊺_n+iRu_n+i) )

. (2.5)

The matricesQ and Q_y hold a relation

Q=C^⊺Q_yC, (2.6)

where C is a matrix with each possible measurement stacked in a way that it pro- duces all the measurable outputs y in the order they are given in 2.5 when right multiplying it with the vector x of the states. This relation can easily be seen by inserting y=Cxinto 2.4.

The minimization in 2.4 is done with respect to choice of control signal at each time step and the choice of measurement on the next time step. The choice of control signals will directly affect the states on the next time step dictated by the process control dynamics, and the measurement at the next time step will affect the quality of information before next control action is chosen.

2.4 Separation of the cost function

Next we show that the cost function 2.4 related to the control problem can be separated into a sum of two parts. One of the parts will consist of taking minimum with respect control signal only and the other part will consist of taking minimum with respect to measurement choice only. This property on an LQG-control problem is called the separation principle. The benefit of the separation is for example that in the use cases in this thesis, we will assume that the plant already has implemented

an optimal controller for the process, but their options for measuring the process outputs are limited. Then we can assume that the part of the cost function that deals with controls has been minimized. Hence we need only to deal with the part consisting of minimizing with respect to the measurement choices. To prove this separation of control and measurement, we need the next theorem about the expected value of a quadratic form:

Theorem 1. Let x: Ω→Rⁿ be a normally distributed random variable x∼N

µ,Σ

and H ∈R^n×n a symmetric matrix. Then it holds that E{x^⊺Hx}=µ^⊺Hµ+ tr(HΣ) Proof. (see Bates 2010, Thm. 6)

Now we can prove the main result. The proof will be done by induction with respect to optimization horizon N. Similar proof can be found in (Meier, Peschon, and Dressler 1967).

Theorem 2. LetV be the cost function of an LQG-control problem as in 2.4. Then for any N ∈N it holds that

V_N

µ_n+,Σ_n+

=µ^⊺_n+Q_Nµ_n++q_N(Σ_n+) for some matrixQ_N and function q_N.

Proof. First, we prove this for N = 1:

V₁

µ_n+,Σ_n+

= min

E{x^⊺_n+1Qx_n+1}+u^⊺_nRu_n

= min

µ^⊺_n+1−Qµ_n+1−+ tr(QΣ_n+1−) +u^⊺_nRu_n

= min

(Aµ_n++Bu_n)^⊺Q(Aµ_n++Bu_n) +u^⊺_nRu_n

+ tr(QΣ_n+1−) min

(u_n+K₁µ_n+)^⊺Φ₁(u_n+K₁µ_n+)

+µ^⊺_n+Q₁µ_n++ tr(QΣ_n+1−).

Above we used Theorem 1 along with the predict step of Kalman filter. The last equation is how we want to present our cost function. There the term inside the minimization will reach zero by choosingu_n wisely, and hence the cost function will be of the claimed form. The gain term K₁ is called the optimal feedback gain, and

it gives the optimal controlu_n =−K₁µ_n+ when we have an estimate of the current state. It is noteworthy, that the term consisting of the estimate variance is taken out of the minimization. This is because no measurement has yet occurred at this time. Normally this would have to be included in the minimization, since in the future the variance will depend on measurements chosen.

In order to show the claim for N = 1, we need to show that the parameters Φ₁ and Q₁ introduced in the last form of V₁ actually exist. Since the minimization term is zero by choice of u_n, and Σ_n+1− is uniquely determined by Σ_n+ through the predict step of Kalman filter and process dynamics, the cost function will be of the claimed form if these parameters exist. By comparing the terms that are left multiplied byu^⊺_n and right multiplied byu_n, we get

Φ₁ =R+B^⊺QB. (2.7)

Furthermore, by comparing the terms that are left multiplied by u^⊺_n and right mul- tiplied by µ_n+, we get

K₁ = Φ⁻₁¹B^⊺QA. (2.8)

Finally, by comparing the terms left multiplied byµ^⊺_n+ and right multiplied by µn+, we get

Q₁ =A^⊺

Q−QB(R+B^⊺QB)⁻¹B^⊺Q

A. (2.9)

We have omitted the transpose signs of Q, since it is assumed to be symmetrical.

Let us now prove the induction step. Assume that N > 2 and that the claim holds forN−1. This means there exists a matrix Q_N−1 and a function q_N−1 ofΣ_n+

such that

V_N−1

µ_n+,Σ_n+

=µ^⊺_n+Q_N−1µ_n++q_N−1(Σ_n+).

We want to show that the claim holds forN. The trick is to separate the first step in the cost function and use the assumption for the rest of the horizon. We shall define a new notation for the minimization, wheremin_n→ means minimization from time stepn to the end of the horizon, and min_n means minimization wrt. just one

time step n. With this notation we get

V_N

µ_n+,Σ_n+

= min

n→ E (_N−1

X

i=0

(x^⊺_n+i+1Qx_n+i+1+u^⊺_n+iRu_n+i) )

= min

n→ E

x^⊺_n+1Qx_n+1+u^⊺_nRu_n

+E (_N−1

X

i=1

(x^⊺_n+i+1Qx_n+i+1+u^⊺_n+iRu_n+i) )!

= min

n E

x^⊺_n+1Qxn+1

+u^⊺_nRun

!

+ min

n min

n+1→E (_N−2

X

i=0

(x^⊺_n+i+2Qxn+i+2+u^⊺_n+i+1Run+i+1) )

= min

n E

x^⊺_n+1Qxn+1

+u^⊺_nRun

!

+ min

n VN−1

µn+1+,Σn+1+

= min

n E

x^⊺_n+1Qx_n+1

+u^⊺_nRu_n+Ey

µ^⊺_n+1+Q_N−1µ_n+1+

+q_N−1(Σ_n+1+)

!

The order of minimization and summation can be interchanged like they are in the third equality, since all the terms are non-negative due to Q and R being positive semi-definite. Also in the third row the terms within the first minimization do not depend on future control or measurement choices, hence the minimization is taken only wrt. current time. The expected value with subscripty in the formula ofV_N−1 in the time stepn+1means that that the meanµ_n+1+ depends on measurement data y_n+1that is not yet obtained, but will be available before the controlu_n+1+is chosen, and that the measurement data y_n+1 is normally distributed. Hence the estimate mean will be normally distributed as well according to Kalman filter update step.

According to Kalman filter, the variance Σ_n+1+ is deterministic, no matter what kind of measurement data is obtained at the time stepn+ 1. For this reason, only the expected value of the terms containing µ_n+1+ must be taken. We have already solved the first expected value in the case of N = 1, but now we must find out the expected value and variance of µ_n+1+ in order to use Theorem 1. From Kalman filter we have

µ_n+1+=µ_n+1−+L_n+1

y_n+1−C_n+1µ_n+1−

and furthermore, since the expected value of the measurement datay_n+1isC_n+1µ_n+1− and the variance of the estimate evolves in a deterministic way according to Kalman filter independent of the measurement data obtained, we have

µ_n+1+ ∼N

µ_n+1−,Σ_n+1+

.

Finally after applying Theorem 1 to the two expected values of quadratic form in

the last form ofV_N, and Kalman filter predict step backwards, we have V_N

µ_n+,Σ_n+

= min

un

(Aµ_n++Bu_n)^⊺Q(Aµ_n++Bu_n) + (Aµ_n++Bu_n)^⊺Q_N−1(Aµ_n+Bu_n) +u^⊺_nRu_n

+ min

Cn+1

tr(QΣ_n+1−) + tr(Q_N−1Σ_n+1+) +q_N−1(Σ_n+1+)

.

(2.10) The terms have been split into those that depend only on control and to those that depend only on measurement. The index of measurement isn+ 1since this is when the first measurement is made. Now the functionq_N is the second minimization of the equation above. We still need to show that there existsQ_N of the wanted form.

We write the terms inside the minimization wrt. control in the form (u_n+K_Nµ_n+)^⊺Φ_N(u_n+K_Nµ_n+) +µ^⊺_n+Q_Nµ_n+,

since if the parameters Φ_N and Q_N exist, the minimum wrt. control will be u_n =

−K_Nµ_n+ and the claim is finally proven. Comparing the terms like we did in the case of N = 1 we obtain

Φ_N =R+B^⊺(Q+Q_N−1)B K_N =Φ⁻_N¹B^⊺(Q+Q_N−1)A Q_N =A^⊺

(Q+Q_N−1)−(Q+Q_N−1)B(R+B^⊺(Q+Q_N−1)B)⁻¹B^⊺(Q+Q_N−1)

A (2.11) Now we can obtain the claimedQ_N andq_N recursively fromQ₁ andq₁, respectively.

These were shown to derive directly from process dynamics and cost function weight- ing matrices. Hence the claim is proven.

When an optimization horizon N is chosen,QN is easily iterated from its recur- sive formula. After this, the calculation of KN is simple and we have our control law for any given state of the process: the control action at any time stepn is

u_n=−K_Nµ_n+, (2.12)

where µ_n+ is the current estimate mean of the state. This follows from the fact that the minimization wrt. control is formulated in such a way that the choice ofu mentioned above minimizes it without exception.

The minimization wrt. measurement is not so simple. We need to know what is

q_N. For example if N = 3, we get q₃(Σ_n+) = min

Cn+1

tr(QΣ_n+1−) + tr(Q₃₋₁Σ_n+1+) +q₂(Σ_n+1+)

!

= min

Cn+1

tr(QΣn+1−) + tr(Q3−1Σn+1+) + min

Cn+2

tr(QΣn+2−) + tr(Q3−2Σn+2+) +q1(Σn+2)

!

= min

{Cn+i}³i=1⁻¹

tr(QΣ_n+1−) + tr(Q₃₋₁Σ_n+1+) + tr(QΣ_n+2−) + tr(Q₃₋₂Σ_n+2+) + tr(QΣ_n+3−)

!

= min

{Cn+i}³⁻¹i=1

tr(QΣ_n+3−) +

3−1

X

i=1

tr(QΣ_n+i−) + tr(Q_3−iΣ_n+i+)

! .

The nested minimizations above can be combined like this since measurements fur- ther in time do not affect anything in the past. Since nothing changes in the recursion with largerN, we conclude that

q_N(Σ_n+) = min

{Cn+i}^Ni=1⁻¹

tr(QΣ_n+N−)+

NX−1 i=1

tr(QΣ_n+i−)+tr(Q_N−iΣ_n+i+)

!

(2.13) The effect of the measurement choicesC_nare seen in the variancesΣ_n+, which evolve according to Kalman filter determined by the chosen measurements and process dynamics.

Equation 2.13 does not have a simple analytical solution for the measurement sequence like in the case of control. The minimization can be solved off-line for givenN and starting variance Σ_n+ by checking each sequence of measurements and finding the one that minimizes the cost. The measurement sequence can be shown to follow are recurring pattern, and the optimal solution can be solved offline rather than recalculating it at each time step based on current information state.

There may exist ways to reduce the computational burden of the problem by pruning some of the sub-optimal sequences. Alternatively, there exists iteration algorithms that only approximate the optimal solution, but are much faster than a comprehensive search checking all possible sequences (see Ross et al. 2008). The question of what is the best way of solving Equation 2.13 is not within the scope of this study, and for the applications we shall use comprehensive search.

3 Test cases

In this chapter we present the processes studied in this thesis, as well as the mea- surement systems in them. Both of the test cases for measurement optimization have a primary measurement that is the focus of the optimization. The second case has and a secondary measurement that is used for example to detect non-ordinary behaviour of the process that may go unnoticed by the primary measurement, and therefore cause the control to act in a sub-optimal way. The froth flotation process and its primary measurement system, the Courier analyzer, will be introduced first.

After this, we will introduce the process of paper machine basis weight control and its measurement systems, the beta radiation based QCS-scanner and optical camera based web inspection system (WIS).

3.1 Froth flotation process

The separation of control and measurement of an LQG-control problem was tested in the case of froth flotation. A simulator called HSC Chemistry 9 was used in this thesis instead of the real process.

3.1.1 HSC-simulator

The devopement and test environment for froth flotation was the HSC Chemistry 9-simulator, instead of an actual flotation process. The HSC-simulator is Outotec’s process simulator used for chemical reaction calculations and simulations. The sim- ulator can be used to run dynamic and real time simulations and has been developed since the 60s. A simulation of an eight-tank froth flotation line from (Kortelainen 2019) was used to collect data corresponding the actual froth flotation process.

3.1.2 Process and control

In the froth flotation process, a slurry of minerals flows through several flotation tanks filled with water. The process diagram is depicted in Figure 3.1. The slurry entering the first tank is called the feed. The feed is the raw material that will be enriched throughout the process, and may vary considerably in its mineral content.

The tanks contain reagents that for example make minerals hydrophobic. Air is blown through the bottom of the flotation tanks, making bubbles. The hydrophobic minerals will then attach within these bubbles and rise to the layer on top of the tank, called the froth layer, and overflow out of the tank. The overflown froth from each tank is collected into a product called concentrate. Not all of the slurry will be collected in any single tank, and the tanks will flow it into the following tanks.

Figure 3.1Diagram of flotation process. Rectangles depict where grade can be measured.

The slurry that remains after all the flotation stages is called tailings. The flotation process may be continued in different froth flotation lines that may contain different reagents to further serve the purpose of the process. For example, in the following phases the concentrate or tailings of the previous flotation line may be used and the unwanted minerals might be made hydrophobic, and the wanted minerals remain throughout all the following flotation tanks. In the scope of this thesis, we will study eight flotation tanks and only single flotation line.

The main variable of the different phases of the slurry (feed, concentrate and tailings) is its grade, which is the mass of the wanted mineral divided by the mass of the total solids in the slurry. High grade means high purity of the wanted mineral and a part of the goal of the process is to maintain a sufficient grade of the end product to meet some given needs. It is also important to recover as much of the wanted mineral as possible from the feed. The ratio of the wanted mineral in the concentrate to the wanted mineral in the feed is called recovery, and like grade, the higher the recovery, the better the performance of the process.

The grades and recovery of the process are controlled using so called froth speeds.

These are the speeds at which the froth layer on top of the flotation tanks moves out of the tank. The amount of froth collected depends on the diameter of the lip of the tank from which the froth flows and the height of the froth layer, but these are kept constant in this test case. In the HSC-simulator, froth speeds are controlled using PID-control with the air flow to the bottom of the tanks as the control signal.

There exists a tradeoff between the controllable outputs. Generally increasing the

froth speeds of the tanks yields a higher recovery, but lower concentrate grade. We define two input signals for the process: FrothSpeed 1, which is the setpoint of the controller of froth speeds of the first four tanks, and FrothSpeed 2, which is the setpoint of the controller of froth speeds of the final four tanks.

The froth flotation process is very complex and has many intricate variables that affect each other in various ways. In this thesis, we will not go deeper into these intricacies, but instead will analyze the process as a simple input-output system.

A deeper analysis of the actual physics and chemistry of the froth flotation process can be found for example in Fuerstenau, Jameson, and Yoon 2007. A linear state space model of the process was built. This model was used to determine the optimal measurement sequence for the process. The modeling, optimization and testing of the froth flotation process will be covered in Chapter 4.

3.1.3 Measurements

The primary measurement system of the froth flotation process is an X-ray fluo- rescence (XRF) based Courier analyzer. The Courier analyzer can make on-line measurements of the grades in different streams of the froth flotation process ex- plained in Section 3.1.2. The Courier analyzer can only sample and measure the grade in one stream at a time in the scope of this study. Measurement time may be changed according to needs. Short measurement time ensures fast reaction to possible changes in the process, but comes at a cost of less accurate measurement.

Longer measurement time can be used to give more accurate measurements, but then the process may change during the measurement and some changes may go unnoticed (Remes, Saloheimo, and Jämsä-Jounela 2007). The focus of this study is to determine which stream to measure with the Courier analyzer at each time instant to optimize the process control, rather than how long the measurement time should be.

3.1.4 Test case objective

The grades of the slurry in different streams are measured using the Courier analyzer.

The problem arises when there are not enough Courier analyzers to measure the grade in all the different streams. In this thesis, where a single froth flotation line is studied, there are 3 different possible measurements: feed grade, concentrate grade and tailings grade. The LQG measurement optimization problem then becomes to determine which measurement to choose at which time instant when the optimal LQG-control is used. The problem can be extended to consider multiple different froth flotation lines, but this is out of the scope of this thesis and will require either much more processing power or better search algorithms to solve for higher

Figure 3.2 A diagram of the paper machine.

optimization horizons.

3.2 Paper machine basis weight control 3.2.1 Process and control

As another application of optimization of the information state, the basis weight (grams per square meter) of a paper machine was considered. A rough diagram of the process can be seen in Figure 3.2. The pulp slurry from the short circulation of a paper machine arrives at the headbox, where it is distributed as a wide (around 4-11m) flow through a slice opening, to be gradually dried to a paper web. This flow travels at high speed (around 10-30m/s) through pressing and final drying and is reeled as a finished paper product. The objective of the basis weight control is to ensure uniform quality paper across the entire web at all times. Cross-directional (CD) variation of basis weight means a spatial deviation from target, persistent in time. Machine directional variation (MD) means deviation in time, uniform across the web width. An example of a CD disturbance is that there is a persistent bump in the middle of the web in relation to the rest of the width of the web, whereas basis weight is oscillating in time throughout the width of the web due to inconsistencies in the incoming pulp slurry is an example of MD disturbance.

The basis weight is controlled in CD using slice opening or dilution actuators, tens of which are placed along the width of the headbox where the web is formed.

Larger slice opening in a point of the CD means more pulp slurry will enter the

web at that point, which in term increases the basis weight locally. Dilution water actuators add more water to a given point in the CD, decreasing the basis weight.

The CD-control of the process will take action on deviations from the average basis weight of the web in a certain location in the CD of the web. Disturbances in CD are spatially local and corrected by the CD-control. When a disturbance in the pulp slurry consistency occurs before the headbox, it will be temporal in nature.

This is to say that these kind of disturbances affect the entire width of the web.

Temporal (or MD) disturbances are corrected using MD control, which aims to keep the average basis weight of the web within the set standards.

3.2.2 Measurements

The primary measurement of the web basis weight is done by a QCS-scanner moving across the width of the web. The scanner measures the basis weight using trans- mittance of beta radiation. The scanner can only measure very small portion of the web, and moves slowly (around 0.2-0.5m/s) compared to the speed of the web in machine direction. For these reasons, only a portion of the information state of the web can be updated at any time instance with the scanner. The area measured by the scanner at a single time instance is typically around 1 centimeter wide and it takes typically up to several tens of seconds until the same position is measured again. Hence the utilization of the scanner measurement across the wide web is a critical question. If moving at a constant speed, the scanner will measure the basis weight in a zig-zag pattern since the entire web is moving in time. The speed of the scanner can be varied to measure certain locations more or less thoroughly, and the scanner can even be stopped to measure variation in time in a certain location.

The measurements of the scanner produce and estimate of the basis weight along the entire width of the web. The estimate is updated using a Kalman filter. Since the measurements made by the scanner are very limited spatially, the uncertainty of the estimates is usually very high. A detailed analysis of the scanner control can be found in (Raunio and Ritala 2018).

The basis weight can also be obtained from the optical web inspection system (WIS). The WIS is a camera-based system that is meant for detecting larger flaws in the web, such as holes or dirt particles, but may be used to measure the basis weight to some extent. It consists of several cameras and an illumination panel of LEDs. The images can be used to measure optical transmittance by inspecting the light scattering from the air-solid interface of the web. Optical transmittance in turn can be used to obtain the basis weight. The problem is that the scattering of the light depends on many other factors than basis weight also. For example moisture, filler content and color of the paper affect the light scattering, and hence the optical transmittance. For this reason the basis weight cannot be controlled only with the

information provided by the WIS and the more accurate scanner measurements are necessary.

The primary information state optimization of the process will come from schedul- ing the scanner according to current information state and to some extent, the infor- mation obtained by the transmittance measurement. More about the process that is the basis weight control of a paper machine can be found in (Ritala 2009). The chapter of the book covers the paper machine process control extensively.

3.2.3 Process model

The process will be analyzed as a linear, almost integrative state space model with disturbance terms (see Raunio and Ritala 2018, p. 78)

x_n+1 =Ax_n+B_cd∆u^CD_n + ∆u^{M D}_n−d

md+D_n+1+ϵ_n+1. (3.1) The matrix A is nearly integrative. The states themselves are a partition of the width of the web at the headbox. This means that each state is a small portion of the web in CD. The finer the partition, the more accurate the model. As a downside, more states means more complex calculations.

The matrixB_cdis the response model for the slice opening actuators that control the basis weight in CD and∆u^CD_n is the change in the CD actuators between time steps n −1 and n. The steady state unit step response of a single slice opening actuator can be seen in Figure 3.3. A slice opening actuator will increase the total flow into the web in its center, but decrease the flow near it. The change in actuator may be instantaneous, or first order transfer function dynamics may be introduced to give the response a rise time.

The term∆u^{M D}_n is similarly the change in MD actuator affecting the entire width of the web between time stepsn−1andn and contains first order mixing dynamics before the headbox. There is also a pure delay ofd_md between the MD actuator and the headbox.

The termsDandϵare load disturbance and process noise respectively. The load disturbance will be added to test different scenarios and can be a mixture of CD disturbances and MD disturbances. The process noise is Gaussian white noise in time, possibly correlated in CD.

The two measurement systems of the process model are the QCS scanner and WIS. Neither of these measurements are made at the headbox, but further down in MD, meaning that the measurements themselves are delayed wrt. states. The delays of the two measurements are not exactly the same, but of the same order of magnitude, in this analysis 55s for optical measurement and 60s for scanner.

The information state is updated according to the scanner data. Estimates are

Figure 3.3 Steady state unit step response of a slice opening actuator in the middle of the web.

built using a delayed Kalman filter that takes into account the delay between the states and the measurement (see Raunio and Ritala 2018, p. 79). The predict step of the delayed Kalman filter is like that of a normal one, except that the delay is taken into account when adding the actuator effects. The delay is taken into account by first applying Kalman filter to the delayed estimates, and then predicting to the present time using the predict step as many times as there are time steps in the delay. The component-wise predict step of the delayed estimate is

µn−d+1−[i] = µn−d+[i] +

acts.

X

j

Bcd[i, j]∆u^CD_n−d[j]

Σn−d+1−= Σn−d++ Σ^(p),

(3.2)

wheredis the delay from the headbox to the scanner. The matrixAis left out since it is very close to identity in this use case and will not have much of an effect on the estimates or their uncertainty. The summation in the first row of 3.2 refers to summing over each of the CD actuators. The component-wise measurement update step of the delayed Kalman filter is

µ_n−d+1+[i] =







1−α_n−d+1[i]

µ_n−d+1−[i] +α_n−d+1[i]y_n+1[i] if i=i^(scan)_n+1

µ_n−d+1−[i] else

Σ_n−d+1+[i] =







Σ⁻_n−¹_d+1−[i] + Σ^(m)_n+1⁻¹ ₋1

if i=i^(scan)_n+1 Σ_n−d+1−[i] else,

(3.3)

where y_n+1 is the new measurement data obtained, i^(scan)_n+1 is the scanner index CD position and

α_n[i] = Σ^(m)_n+1⁻¹ Σ⁻¹_n−[i] + Σ^(m)_n+1⁻¹

for all i. (3.4)

The delayed Kalman filter estimate and its uncertainty are only updated in the CD positions where the scanner is. In other areas, the measurement update step will do nothing. The update step for the estimate becomes an exponential filter with the weighting factor α_n depending on the uncertainty of the estimate at the specific time instant. If the uncertainty is high (ie. it has been a long time since the last measurement), then more weight will be placed on the new measurement data obtained rather than on the old estimate. The estimate and its uncertainty are then predicted to the current time using Kalman filter predict step recursively dtimes:

µ_n+1+[i] =µ_n−d+1+[i] +

acts.

X

j

B_cd[i, j]

u^CD_n [j]−u^CD_n−d[j]

Σ_n+1+= Σ_n−d+1++dΣ^(p)

(3.5)

The first equality follows from the fact that the matrix A is close to identity and hence left out, and that the control terms u^CD_n between the time steps that are left are cancelled out in the recursion.

The optimal CD control law is such that it attempts to fix any CD deviations from the mean of the estimate in steady state according to given cost function

V = min

∆u^CD_n

¯

x^⊺_n+1Q¯x_n+1

, (3.6)

wherex¯is the component-wise deviation of the state vector from its target profile.

The state weighting matrix Q of the quadratic form is chosen to be identity in this test case. At each time step, optimal change in control ∆u^CD_n is a such that

it attempts to counteract the deviations from the mean using the response model matrixB_cd of the slice opening actuators. Assuming the disturbances and noise are zero, we get from equations 2.12 and 2.8 that

∆u^CD_n =− B^⊺B₋1

B^⊺µ¯_n, (3.7)

where µ¯_n is the component-wise deviation of the current estimate from its target profile. If the change in actuators has a rise time rather than being instant, only a portion of the control should be implemented in order to avoid stability issues.

Where CD control takes care of any deviations from the target CD profile, the point of the MD control is to regulate the CD mean to a given setpoint. When calculating the mean basis weight, the states can be given weights proportional to the inverse of their current uncertainty. The weighted mean will then depend less on the state estimates whose uncertainty is high and more on those whose uncertainty is low. The control law in this study is that of a PI-controller:

u_n =−K(e_n+ I_n T_i) I_n+1 =I_n+e_n,

(3.8)

whereK is the controller gain, e is the deviation of the mean basis weight from its setpoint, I is the integrator state value and T_i is the integration time.

3.2.4 Process simulator

The tests were run of the process simulator in Matlab. The simulator consists mainly of definitions, storages, and the actual iterative simulator. In the definitions, the user defines parameters for the process, such as web width and its partitioning into discrete states, number of CD actuators, scanner speed, web speed, disturbances, level of control, etc. Storages exist for storing simulated data such as real-life sim- ulated data, estimates and different profiles. The states of the process are small cross-directional sections of the web at the headbox. In the test cases, a 5 meter wide web will be partitioned into 0.5 centimeter wide sections and this means the process will have 1001 states. For this reason, the storages will need to be buffered and the final data downsampled so that the simulator can be run for longer periods of time without running out of memory.

The actual simulation is an iterative loop with respect to simulation time. Mea- surement data along the width of the web is simulated and stored according to defined parameters and process model. Disturbances and noise are added to the simulated data, as well as control action. Web shrinkage, defined as shrinking of

the basis weight in the middle of the web, may also be added to the real-life data.

Optical transmittance measurements are then simulated across the entire width of the web and stored. The scanner motion is then simulated, as well as measurements made by the scanner in its CD position. The scanner may move in one time step across multiple states, and the measurements taken will then be averaged and stored in the storage of the states which the scanner had passed over during that time step.

Estimates of the basis weight along the width of the web are stored according to the scanner data and updated using the delayed Kalman filter. At the end of each scan period (ie. when the scanner has travelled from one edge of the web to the other end.) several profiles of the web will be generated and stored. These profiles are made out of for example raw scanner measurement data collected during the scan period and will display that as a CD profile. At the end of each scan period, the scanner will leave the web area for a period of time.

Finally, the simulator will calculate control action for CD and MD according to the estimates generated from the scanner data according to the control laws introduced in Section 3.2.3.

3.2.5 Test case objective

One of the key problems in the basis weight control is to distinguish between CD variations and MD variations and so called residual variation, which does not fall in either of the first two categories. There are cases where pure MD variations might appear to be CD variations. In this case, the CD-control will act on a deviation in the web that is not spatial in nature, which will cause complications. A part of this study is to figure out ways to detect MD variations that are perceived as CD variations.

The goal is to find solutions for some of the problematic disturbances in real life basis weight control process as well as to optimize the quality of the information state in different scenarios and to test these with the simulator.