On Predicting Milk Yield and Detection of Ill Cows

(1)

On Predicting Milk Yield and Detection of Ill Cows

N A N C Y G R E N N S T A M

Master's Degree Project

Stockholm, Sweden 2005

(2)

Abstract

The fully automated milking system VMS has different functions which com- plements the actual milking of cows. This master thesis presents a method to improve the calculation of milk yield in dairy cows for the VMS. This report also investigates if it is possible to improve the algorithm for finding cows with mastitis (udder inflammation). The correctness of the prediction of milk yield is important for a couple of actions in the VMS. For example, valuable time can be saved if teatcups are attached first to high yielding teats. Only cows with an attained minimum level of predicted yield should be allowed to enter the VMS and get milked. Milking has traditionally been an event to monitor the condition of the cows. Therefore methods that determine the condition are demanded for any automatic milking systems. Mastitis is a costly illness and a working test for ill cows should be implemented in the VMS in order to know which cows that are ill.

The goal of this thesis work is to develop two new algorithms for the VMS.

First, an improved algorithm for the prediction of secretion rate is presented.

The improved algorithm uses a Kalman-filter to update the secretion-rate. The improved method has a lower total prediction in most cases. The Kalman-filter was tested and developed for five farms and was verified on one farm.

Second, this report investigates if a cusum test can be used to detect ill cows.

The method turns out to be slightly better than the current algorithm. A test for cows which are milked on three or two teats is evaluated. In this test the number of milkings with high conductivity and low secretion rate are weighted together. This algorithm is slightly better than the current algorithm used for detection of ill cows.

(3)

Legend

Kick-oﬀ - When a cow kicks the teatcup oﬀ its teat

Lactation - The period after birth of a calf when the cow yields milk Secretion rate - The pace with which a cow produces milk

VMS - Voluntary Milking System, an automated milking system from DeLaval Mastitis - Udder inﬂammation

Mammary gland - Milk producing gland

(4)

Introduction

1.1 Background

1.1.1 The master thesis

This master thesis project consists of two parts and has been conducted at DeLaval International AB. DeLaval department of VMS (Voluntary Milking System) has a robot which automatically milks cows. For VMS two algorithms associated with the robot have been developed. The ﬁrst part of the project was to ﬁnd an algorithm to calculate Expected Milkyield and the second was to work on a method to detect if the milked cows were ill based on the information gathered in the VMS.

1.1.2 Milking dairy cows (history)

Milk is, and has been for a very long time, an important sustenance for humans worldwide. The production of milk though, has changed quite radically over the last decades. At least in industrial countries milk is produced and sold as any other commodity today, thus increasing the demands on the cows as well as the technique used for milking the cows.

The average milk production 30 years ago of a cow was 4000 kg per year.

Genetic progress as well as increased knowledge of heard management has made the average cow of today yield between 7000 and 12000kg milk per year [1].

Automatic milking is a conception used for the technology of milking the cows with a machine in contrast to milking by hand. Still, when using the traditional automatic milking system one has to attach the teatcups to the cows’ teats manually, and with increasing heard sizes this is a labor that get more and more physically demanding. But lately, alternatives to these systems have come available for dairy farmers.

1.1.3 Voluntary Milking System by DeLaval

Since 2000 DeLaval has manufactured a milking system in which also the attaching of the milking machines’ teatcups is fully automated. In traditional automatic milking the operator decides when the cow is milked, i.e. when the

(7)

basis. The cow decides for itself when it wants to be milked. DeLaval calls their system Voluntary Milking System, VMS.

The VMS includes a milking station, MS, where the cow is milked. Unlike traditional milking, every teat is milked and monitored individually in VMS, and since all cows that enter the MS are identiﬁed by the system the milking procedure can be tailored for each cow and teat.

1.1.4 Expected milk yield

The VMS renders the possibility to make the milking procedure more eﬃcient.

By predicting the amount of milk the cow will yield, the milking can be adapted for that amount of milk.

Several actions can be based upon the expected milk yield.

1. A cow which predicted yield is over a certain level should be allowed to get milked.

2. By attaching the teatcup on the highest expected yielding teat ﬁrst, valuable time can be saved in the MS.

3. If a cow kicks the teatcup oﬀ the teat, there is no point in reattaching it if the cow has very little more milk to give.

4. If a cow yields less milk than expected it could be a sign of illness and should be noted.

These and more applications rely on the correctness of the expected milk yield calculation.

An estimate of the amount of milk a cow will yield can be derived from the knowledge of prior milk yield amounts and the time since the cow last was milked.

1.1.5 Ill cows

Milking has traditionally been a time to monitor the cows’ conditions and health.

If it is possible to draw conclusions on the cows’ health from data that could be collected in the VMS, then such conclusions should be drawn. Illnesses such as mastitis need to be found and treated to ensure an eﬃcient milk production.

1.2 Thesis purpose

This thesis is divided in two main parts. The ﬁrst is partly done together with Nicklas Lundin and considers the calculation of expected milkyield in VMS. The goal is to present a proposal for an algorithm that can be implemented in the VMS. The proposed method should use the information of the cows changed milking characteristics throughout the lactation. The performance of the milk yield predictor should be evaluated against the one currently used in VMS by comparing the actual milk yield with the predicted. The second part of this thesis deals with the detection of ill cows. The farmer needs information about which cows that are ill in order to treat them. The detection should be based on information gathered in the VMS, primary conductivity and secretion rate.

(8)

Part I

Expected Milk Yield

(9)

Chapter 2

Theory

2.1 The lactating cow

2.1.1 Milk yield between milkings

The velocity with which the cow produces milk is called secretion rate, measured in gram per hour (g/h). The secretion rate x is not a linear function, i.e the total yield is not equal to x∆t, where ∆t is the time since last milking. The yield is for example saturated after some time with a time which is individual for each cow. This is because of the limited size of the cows’ udder. It cannot carry an unlimited amount milk.

The literature on lactating cows is somewhat vague concerning how the secretion rate declines. After 10 hours the secretion rate start to decline and after 20-35 hours it has stopped [3]. The milk contains a chemical inhibitor which suppresses the milk secretion when the amount if milk is increasing. If the milking interval is over 16h this inhibitor can aﬀect the secretion rate in the next milking as well, at least up to 8h after the milking.

2.1.2 Milk yield throughout a lactation

The ﬁrst time after the birth of the calf the cow undergoes rapid changes in milk secretion rate. When the calf is born the secretion rate of the cow increases for approximately four weeks until it hits its peak [1]. How the secretion rate changes after that is very individual, but it is of course favorable to have the cow stay at this peak level for as long as possible. The time after this peak usually incorporates a slow decline in secretion rate until the cow goes dry, i.e.

the end of her lactation. A typical lactation is shown in Figure 2.1.

2.2 Method currently used in VMS

The method currently used for calculating the expected milk yield uses a auto regressive ﬁlter when updating the secretion rate, and it uses a time transform equation to compensate for non-linearities in the secretion rate. To get the expected milk yield, the secretion rate is multiplied with the time since last milking, ∆t.

(10)

3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 x 10⁴ 0

100 200 300 400 500 600

Yield per hour [g/h]

time [h]

Figure 2.1: The secretion rate for one teat over one lactation is shown.

2.2.1 Time transform

To compensate for the saturation in milk secretion described in section 2.1.1, this method uses a time transform. When ∆t > 16h the time between the adjacent milkings, ∆t, is transformed as if the secretion rate was linear.

If 16h < ∆t < 35h the transformed time ∆t_transis calculated as

∆t_trans= a₃(∆t− a0)³+ a₂(∆t− a0)²+ a₁(∆t− a0) + a₀ (2.1) For ∆t < 16, ∆t_trans= ∆t.

2.2.2 Auto regressive ﬁlter

The secretion rate x(t) is estimated as ˆ

x(t + 1) = ∆t

T_constx(t) + ˆx(t)(1− ∆t

T_const) (2.2)

where Tconst is a constant.

This method does not take into consideration where the cow is in its lactation since T_constdoes not change.

2.2.3 Expected yield

Finally the expected milk yield ˆy(t+1) is calculated as ˆy(t+1) = ˆx(t+1)∆ttrans.

(11)

2.3 The proposed algorithm

The basic idea is to predict the milk yield from the secretion rate and the time since the last milking. A Kalman ﬁlter can be used to update and predict the secretion rate, even if the secretion rate is non-linear.

2.3.1 Kalman Filter

The Kalman ﬁlter estimates the state of a process in a way that minimizes the mean squared error. It can be used to estimate past, current and future states.

In this application it will be used to predict future values of milkyield.

General Model

Consider the system in (2.3). The state x(t) can not be measured. The future values of x(t) are aﬀected with process noise v(t). The output is y(t) and e(t) models the measurement noise.

x(t + 1) = F x(t) + Gu(t) + v(t)

y(t) = Hx(t) + e(t) (2.3)

Kalman ﬁlter

By using the Kalman ﬁlter [2] the prediction of x(t) will be as follows ˆ

x(t + 1|t) = F ˆx(t|t − 1) + Gu(t) + K(t)[y(t) − H ˆx(t|t − 1)] ,

The feedback gain is calculated as follows

K(t) = [F P (t|t − 1)H^∗+ R₁₂]Λ⁻¹(t) ,

where

Λ(t) = HP (t|t − 1)H^∗+ R₂ ,

and

P (t + 1|t) = F P (t|t − 1)F^∗+ R₁− K(t)[R21+ HP (t|t − 1)F^∗] . The constants R₁, R₂, R₁₂ and R₂₁ are

E

e(t) v(t)

_∗

=

R₁ R₁₂ R₂₁ R₂

.

Model for lactating cow

In our system x(t) is the secretion rate (see 2.1.1). As the secretion rate can not be directly measured, y(t) will be the calculated yield per hour.

We have a scalar system.

x, F, G∈ R.

(12)

In the model of the milk secretion process we have G = 0 and H = 1. In the model F = 1 for the ﬁrst milkings since the milk yield is increasing in the beginning.

The noises v(t) and e(t) are modelled as uncorrelated white noise sequences, N(0, R₁) and N(0, R₂). This give R₁₂ = R₂₁ = 0. The process noise range from excluded model parameters to individual characteristics for each cow. The measurement noise is of course measurement errors in the milk meters, but it also includes the time transform error. Since values have to be assigned to R₁ and R₂, they can be seen as design parameters for the Kalman predictor.

The model for the lactating cow will be as follows:

x(t + 1) = F (t)x(t) + v(t)

y(t) = x(t) + e(t) (2.4)

The model (2.4) updates the Kalman predictor gain K(t) and predicts the se- cretion rate ˆx(t) , with initial estimates of ˆx(t|t − 1) = x₀ and P (t|t − 1) = P₀, as

K(t) = F (t)P (t)

P (t) + R₂ (2.5)

ˆ

x(t) = ˆx(t− 1) + K(t)[y(t − 1) − ˆx(t − 1)] (2.6) P (t + 1) = F²(t)P (t) + R₁− K(t)P (t)F (t) (2.7)

Equations (2.5) to (2.7) are then repeated for every milking.

In order to ﬁnd the stationary values for P and K we let P (t + 1) = P (t) and F (t) = 1 as t→ ∞. From (2.5) one gets the constant value for K:

K = P

P + R₂ , (2.8)

and from (2.7) one gets

P = P + R₁− KP , (2.9)

which together gives the positive solution for P

P = R₁ 2 (1 +

1 + 4R₂

R1) . (2.10)

Inserting this equation in (2.8), one gets the stationary Kalman gain

K = 1 +

1 + 4^R_R²

1

1 +

1 + 4^R_R²

1 + 2^R_R²

1

(2.11)

(13)

2.3.2 Non-linear secretion rate

The non-linearity in secretion rate due to saturation in the cows udder, brings a time variance into the state-space model (2.4). But by using a time transformation like the one described in (2.1) the state space model can be used as a time invariant system.

2.3.3 Calculating the actual secretion rate

Since the secretion rate can not be measured directly it is derived from the milk yield. The relation between milk yield, z(t), and secretion rate, x(t), is z(t) =

∆tx(t), but since the trajectory of ∆t is described in the time transformation, Section 2.1, the secretion rate can be calculated as x(t) = _∆t^z^(t)

trans.

2.3.4 Prerequisites

Since it is not possible (for this thesis work) to run tests on an actual running VMS, data from past milkings done on Milking Stations on a number of farms since the launch of VMS in 2000, is used in the development and validation of this thesis work.

In VMS every cow has a file in which data from her milkings are logged. It is possible to follow a cow through her lactations via these files, presupposed that it has been milked in the MS and not elsewhere. These files are called, because of their suffix, res-files. In the res-files the time of milking, actual milk yield, milk yield predicted by the current algorithm and other logged data can be found. It is important to note that the current milk yield predictor may affect the actual milk yield. The predicted milk yield can only affect the actual milk yield when a kick-off occurs. A kick-off is when the cow kicks the teat cup off her teat. If a kick-off occurs when the attained milk yield is over a certain percentage of the predicted milk yield, the MS will not reattach the teatcup thus letting the predicted milk yield affect the actual milk yield. The MS detaches the teatcup when the speed of the milk from the teat is below a defined level and not when the predicted milk yield is attained.

2.3.5 Methods of comparison

In order to evaluate and compare the expected yield calculation algorithms, methods to do so are necessary. Two complementing measurements are used in the development of these algorithms.

To compare two methods the diﬀerence between actual yield and predicted yield, the prediction error, for all milkings is calculated. The ﬁrst method is to simply compare the sum of the prediction errors for all milkings in one lactation in the current algorithm with the proposed algorithm. Hence

total diﬀerence =

number of milkings

k=1

|actual yieldk−predicted yieldk| (2.12)

The second method is to compare how the distribution of the relative errors diﬀer e.g compare how often the two algorithms have relative errors in the different intervals.

(14)

f ≤ 5%, 5% < f ≤ 10%, 10% < f ≤ 25%, 25% < f ≤ 50%, 50% < f where f = |actual yield−predicted yield|

actual yield

With these two methods it is possible to get a basic value for comparison between algorithms and by studying the distribution, one can come to more detailed conclusions.

(15)

Chapter 3

Analysis and Results

3.1 Time transform

Alternatives to the time transform currently used in the VMS, which is described in Section 2.1, have been evaluated. One of them is described below: If 10h <

∆t < 35h the transformed time ∆t_transis calculated as

∆ttrans= a₂(∆t− a₀)²+ a₁(∆t− a₀) + a₀ (3.1) where a₀= 10, a₁ = 1 and a₂ =−0.02. For ∆t < 10, ∆ttrans = ∆t. A plot of this transform is shown in Figure 3.1

0 5 10 15 20 25 30 35

0 5 10 15 20 25

Real time (h)

Transformed time (h)

Figure 3.1: Alternative time transform of ∆t.

This transform follows the theory [3] that the secretion rate starts to decline after 10h. But, as mentioned before, it is hard to ﬁnd an unambiguous descrip- tion of how the secretion rate declines. It is therefore likely to assume that this

(16)

is individual for each cow. It is also likely to assume that it is dependent on the initial secretion rate - the udder will be full sooner if the secretion rate is higher.

The tests done with diﬀerent time transforms indicate that so is the case. No transform giving a better result on all lactations could be found. One transform could give better results on cow A and worse on cow B, when another transform made the other way around.

The important conclusion is that a more detailed time transform is hard to get. It would have to be adapted to the individual cow. The main function of the current time transform is to model the saturation really coarsly.

3.2 Updating the actual secretion rate

Some conditions to the method used to calculate the actual secretion rate, see Section 2.3.3, are added. These conditions include, if the time between two adjacent milking was less than a set time the last secretion rate should be weighted in the calculation as well, since most milkings have a larger ∆t. If two milkings are logged almost immediately after each other it is probably the result of a kick-oﬀ or some other failed milking and the secretion rate cannot be updated in the usual way without getting a very high value. Also if the milk yield is very low there can also be assumed that some abnormalities have occurred.

When ∆t < 6h the actual secretion rate is calculated as x(t) = ^z^k_∆t^(t)+z^k−1^(t)

k+∆tk−1 , i.e the previous milking is accounted for in the calculation.

When ∆tk+ ∆tk−1< 0.1h the secretion rate is not updated at all.

When the milk yield is below 200g the secretion rate is not updated at all.

If the secretion rate is higher than 1000g/h then the secretion rate is not updated.

If the predicted secretion rate deviate more than a set value from the actual value one might not wish to continue to update the secretion rate in the same way. So if |x(t) − ˆx(t)| > θ(x, t) the update has to be altered. θ(x, t) is chosen to be dependent of the standard deviation σ, as θ(x, t) = τ σ(t) and τ = 4.

This method for calculating the tolerance is adapted from the method currently used in the VMS. For the ﬁrst 5 milkings it is initiated as:

σ = 1 n− 1

ⁿ

k=1

(x_k(t)− ˆxk(t))² , (3.2)

where

n = 1, 2, 3, 4 or 5 . (3.3)

Later the standard deviation is updated as a moving average ﬁlter:

σ(t) =

λ(x(t− 1) − ˆx(t − 1))²+ (1− λ)σ²(t− 1) (3.4)

(17)

If the predicted secretion rate fault is greater than the tolerance, then the secretion rate is updated with a smaller Kalman gain, K_new, in feedback. Con- sequently, if|x(t) − ˆx(t)| > θ(x, t) then Knew(t) = 0.45K(t) for that milking.

3.3 Evaluating the Kalman ﬁlter

The performance of the developed algorithm is compared to the performance of the algorithm currently used in VMS.

3.3.1 Parameters in the Kalman ﬁlter

The parameters for the proposed algorithm are presented here. The theoretical background is presented in Section 2.3. Since it is known that the milk yield is increasing in the beginning of the lactation, F should be larger than one in the beginning. But because the behavior of the system in the beginning is more noisy than otherwise F is chosen as only slightly larger than 1. The average number of milkings for which the secretion rate is increasing is about 33. So for the ﬁrst 33 milkings F (t) = 1.0005− 0.000015t. After that F is set equal to 1.

After a number of iterations, or milkings in this case, the Kalman gain settles to a constant K, see Figure 3.2. This constant is chosen small so the perturbations are reduced but great enough to follow the changes in secretion rate throughout the lactation.

The stationary Kalman gain is optimized from simulations done on cows more than six weeks into their lactations. When using a P-controller (equiva- lent to having a Kalman predictor with a stationary gain) it was found that, depending on farm and cows, a feedback gain of 0.2 < K < 0.4 gives the smallest total prediction error. A low K lowers the noise sensitivity and, from equation 2.11, R₁= 1 and R₂= 8 gives a stationary Kalman gain of approximately 0.3

It is hard to predict the secretion rate in the ﬁrst milkings, so K is chosen close to 1 the ﬁrst milkings. This corresponds to P₀≈ 1000.

With

P₀= 1000, R₁= 1, R₂= 8

the Kalman gain K starts at approximately 0.99 and has settled to 0.3 after 12 milkings.

If a cow re-enters the VMS after a time of absence and is not in the beginning of her lactation, it is preferable to have another parameter setting. It would probably not have the same secretion rate as it had when it was taken out of the system. A high K value would allow the predictor to quickly adapt to the current secretion rate, but keeping it for 12 milkings would lead to a very noise sensitive predictor.

3.3.2 Comparing with current algorithm

The algorithm was tested on cows from ﬁve farms. A sixth farm was used for validating the algorithm.

The ﬁles containing data on the milkings for each cow were divided into separate blocks. Each block is meant to represent a cow’s lactation. Since there is no way of knowing (from the ﬁle) where the cow is in her lactation, it is

(18)

0 10 20 30 40 50 60 70 0.2

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

number of milkings

K

Figure 3.2: Feedback gain of secretion rate.

assumed that if it has not been milked for 900h it has given birth to a calf and has started a new lactation.

Lactations which had fewer than 600 milkings were not used in the developing. Also lactations were the cow had been indicated ill were removed. From the farm used for validation only 16 lactations with fewer than 20 milkings were removed.

When using the same time transform as the method currently used in the system, see section 2.1, it is possible to compare the Kalman ﬁlter to the current method.

Table 3.1 shows how many of the lactations where the current algorithm and Kalman algorithm has the smallest total prediction error. The equation for the calculation of prediction error can be found in (2.12).

Farm number Current Kalman Number of lactations

1 7.8% 92.2% 328

2 13.2% 86.8% 552

3 6.7% 93.3% 136

4 32.4% 67.4% 236

6 24.2% 75.8% 496

Table 3.1: Distribution of which method resulted in the lowest total prediction error.

In Figures 3.3-3.5 histograms are used to show the ratio of the total errors of

(19)

proposed algorithm has a lower total error than the current. In Figure 3.3 the ratios of the farms 1-4 together are shown.

−1000 −50 0 50 100

20 40 60 80 100 120

Number of lactations

Ratio between the two methods in percent

Figure 3.3: Distribution of ratios of the prediction errors for the current and proposed algorithm at farms 1-4

In all farms the proposed algorithm had the smallest prediction error in most cases. Farm 4 however is the farm with the lowest percentage of lactations were the proposed algorithm has the lowest prediction error. The ratios from farm 4 are shown in Figure 3.4.

In Figure 3.5 the ratios from farm 6 (the farm used for validation) are shown.

Since farm 4 has quite poor performance compared to the other farms (see Table 3.1), the distribution of the size of the prediction errors on farm 4 is evaluated. The 10 worst (left in ﬁgure 3.4) lactations has a distribution as shown in table 3.2, and the 10 best (right in ﬁgure 3.4) lactations has a distribution as shown in Table 3.3.

f ≤ 5% 5% < f ≤ 10% 10% < f ≤ 25% 25% < f ≤ 50% 50% < f

Moving 47.1% 28.4% 20.3% 2.6% 1.6%

Kalman 47.7% 28.4% 17.9% 2.6% 3.3%

Table 3.2: Distribution of the 10 worst lactations on farm 4

From these tables the conclusion can be drawn that the Kalman ﬁlter has almost the same distribution in the worst case as in the best. This mean that it is the performance of the Moving average ﬁlter that has changed the most. It is also shown that there are some milkings (3.0% in table 3.2) that incorporate

(20)

−400 −20 0 20 40 60 5

10 15 20 25

Figure 3.4: Distribution of ratios of the prediction errors for the current and proposed algorithm at farm 4

f ≤ 5% 5% < f ≤ 10% 10% < f ≤ 25% 25% < f ≤ 50% 50% < f

Moving 37.5% 24.8 % 26.8% 7.2% 3.8%

Kalman 44.3% 28.6% 20.5% 3.7% 2.9%

Table 3.3: Distribution of the 10 best lactations on farm 4

great prediction errors, adding up to the total prediction error. So the Kalman ﬁlter has many milkings with good prediction accuracy but a few with very poor.

(21)

−1000 −50 0 50 100 5

10 15 20 25 30 35 40

Figure 3.5: Distribution of ratios of the prediction errors for the current and proposed algorithm at farm 6

(22)

Chapter 4

Conclusions on expected milk yield

For the milkings to be as eﬃcient as possible in the voluntary milking system the prediction of expected milkyield needs to be made. This thesis has presented an algorithm which predicts the milkyield with higher accuracy than the method currently implemented in the robot.

Cows usually have four teats which are milked. The teats can be seen as in- dependent systems. After a calf is born the cow starts to lactate milk and continues to do so until it goes dry. The milkyield normally increases rapidly in the beginning and declines slowly after the peak is reached. When the milkyield is too low the cow will not be milked any more until the next calf is born. Since the cows decides for themselves when they want to be milked the time intervals between the milkings are not constant. A cow with an empty udder produces milk with a linear secretion rate. The secretion rate will decline as the udder starts to get full. The VMS saves information about when the cows were milked and how much much milk that was gained from each teat. Information on when a cow got a calf is not available.

A Kalman ﬁlter was used to predict the milkyield. The predicted secretion rate is calculated as the previous predicted secretion rate in feedback with a part of the prediction error. The prediction error was the diﬀerence between the previous milk yield and the previous secretion rate. How much of the prediction error that needs to be added to the predicted secretion rate is decided by the Kalman gain. The unlinear part of the secretion rate is compensated with a time transform.

The proposed algorithm for expected yield was in most cases better than the current one. The results vary however at different farms, something which could be understood when taking into account that the different farms operate under different conditions such as different breeds. The farm used for validation shows a result that is worse than most of the other farms’. This could be explained by the fact that this farm did not have as many lactations with insufficient data removed as the other, and that the data at this farm was not used for developing the algorithm. The overall result when validating the algorithm is as good as would be expected. Although the proposed algorithm is not very complicated

(23)

Part II

Ill Cows

(24)

Chapter 5

Theory

The second part of this master thesis deals with developing a new algorithm for the VMS for ﬁnding cows with udder inﬂammation (mastitis).

Mastitis is an inflammation of the mammary gland. Mastitis is a very costly disease of the dairy industry and needs to be controlled. Ill cows should be detected in order to get treatment so the milk production is continually kept at a high level (as mastitis might decrease total milk production) and contaminations to other cows and udders are prevented. The inflammation is the response to infections and irritations in the mammary gland of the cow [3]. The main characteristic of mastitis is the inflammation itself which can vary in severity, depending on both the infectious agent and the cow.

Subclinical mastitis is not detectable by visual inspection as the milk and udder seem normal, whereas clinical mastitis may result in swollen udder and milk clotting as well as in rare cases even death.

As it is of great importance to find cows with mastitis it would be desirable with an improved algorithm. An algorithm for detection of ill cows should have a high sensitivity, which means how many of the ill cows are pointed out as ill. One other variable used for evaluating algorithms is specificity which in detection of ill cows is interpreted as how many of the healthy cows that are considered healthy by the algorithm. The goal for the algorithm is to find 8 of 10 ill cows and not more than 1 of 10 healthy cows should be pointed out as ill.

That means that the specificity should be 90% (of all the healthy cows should the algorithm find at least 90% healthy) and the sensitivity 80% (it should find 80% of all ill cows).

5.1 Mastitis

Mastitis has two major eﬀects on milk; milk reduction and change of milk composition.

Milk from a cow with mastitis has another composition because of the in- ﬂammation. The lactose levels decline which demands ionic changes in the milk. This results in a rising of the concentration of sodium and chloride ions and falling of potassium concentration [3]. This causes the milk to have a higher electrical conductivity. In Figure 5.1 a cow with a high conductivity on the teat

(25)

6.25 6.255 6.26 6.265 6.27 6.275 6.28 6.285 6.29 x 10⁴ 4000

4500 5000 5500 6000

Time [h]

Conductivity [microS/cm]

LF RF LR RR

Figure 5.1: This ﬁgure shows conductivity. This cow had when tested (the solid vertical black line) mastitis on teat left front (LF). The conductivity is on this teat higher than on the other. The teats are named based on the locations of the teats. R stands for rear, F for front, L for left and R for right.

(26)

6.22 6.24 6.26 6.28 6.3 6.32 x 10⁴ 3500

4000 4500 5000

Time [h]

LF RF LR RR

6.23 6.24 6.25 6.26 6.27 6.28 6.29 6.3 6.31 6.32 x 10⁴ 200

250 300 350

Time [h]

Secretion rate [g/h]

Figure 5.2: This cow had when tested (the solid vertical line) mastitis on teat right rear (RR). The upper graph shows conductivity and the lower secretion rate. This cow does not have a high conductivity or a low secretion rate on the teat with mastitis. The teats are named based on the locations of the teats. R stands for rear, F for front, L for left and R for right.

(27)

highest conductivity and the conductivity on this teat is signiﬁcantly higher than the other teats’. But as shown in Figure 5.2 not all ill cows have a high conductivity in the milk from the ill teat and it might not even have a reduced milk production on the ill teat. The conductivity for this cow is not even the highest one, as well as the secretion rate is not even the lowest one for the ill teat. For such a cow that has neither high conductivity nor low secretion rate the detection of the ill teat will be impossible.

5.1.1 How to ﬁnd cows with mastitis

Mastitis can be detected either through the drop in milk yield or through the change of milk composition. The literature states that during a test the relative drop in milk yield at the teats during the infection period was 15.3%± 2.5%

[3]. The average drop in milk yield does not imply that all cows have a drop in milkyield. Also a drop in milk yield does not necessary mean that the cow is ill.

The mastitis should be detected through information gathered in the VMS.

The VMS measures milk yield for other purposes as well and conductivity of milk is also easy to measure. Therefore the idea to ﬁnd ill cows has been to see if a teat has a high conductivity value or/and yields less than otherwise or compared to the other teats.

5.2 The current algorithm used in VMS

The current algorithm for ﬁnding ill cows compares the conductivity value of the milk at each milking. A high conductivity value might indicate that the cow has an udder inﬂammation. But the conductivity value is also dependent on other things for example how much and what the cow has eaten and the amount of water it has drunk. Therefore the current algorithm compares the conductivity value of the teats with each other.

5.2.1 The counter

The current algorithm has a counter for each teat which can take values between 0 and 10. A value over zero indicates that the algorithm has considered this teat ill. If the teat is thought not to be ill the counter goes down one step. The farmer then gets a list on the computer with the cows which are thought to be ill and the value of the counter is also listed. If the value of the counter is rising the value will be colored green. But if the counter is falling the color will be yellow.

5.2.2 The detection

The current algorithm detects mastitis in two ways. They are however both based on conductivity. If a teat which milk has a conductivity value over 15%

higher than the average of the two teats with the lowest conductivity value the cow is ill. If a cow has a conductivity value over 7000µS/cm the teat is considered ill. This is noted and the counter increases by one.

(28)

5.2.3 Results of the current algorithm

In Tables 5.1 - 5.6, it is shown how many of the ill cows the current algorithm found ill in an interval of ±24h. Cows from different farms have been tested at different times. If a cow in the tables is noted ill, it means that the cows were at least one time on the list of the ill cows in the interval. Tables 5.1, 5.2, 5.3 are results from the same farm. The tests are done at different times but not done at the same interval. Tables 5.4 and 5.5 are both the same farm and has fewer cows. One further note on the testing needs to be made. It is not known the exact time of the testings only the date of the test. The tests are assumed to be conducted at midnight. That means that the outcome of the current algorithm is evaluated as how many ill cows the algorithm found the same day as the test was done and the whole day before. Cows from different farms might be of different breeds and be under a different feeding scheme. How and what the cows are fed is not known and the same holds for the breeds. It is however known that certain breeds yield milk with different composition. For example yield Jersey cows less milk but the milk has a higher percentage fat.

This may result in great deviances of any test based on milk composition and yield. A test for indication of mastitis should work for any type of farming and for diﬀerent stocks.

At the ﬁrst test the current algorithm found 30 cows ill . One cow was thought to be ill on two teats. Six of these 30 cows were actually ill and 5 more ill cows were not found by the algorithm, see Table 5.1.

Reality 99 Healthy Cows 11 Ill Cows Outcome 75 Healthy 24 Ill 6 Ill 5 Healthy

Sensitivity 54.5%

Speciﬁcity 75.8%

Table 5.1: The distribution of healthy and ill cows in reality and the outcome of the algorithm for ﬁnding ill cows at test 1.

The second test Table 5.2 is done at the same farm and both sensitivity and speciﬁcity are a little lower. More cows are ill at this test than at the previous (Table 5.1).

Table 5.2: The distribution of healthy and ill cows in reality and the outcome of the algorithm for ﬁnding ill cows at test 2.

Test 3 is also done at the same farm. This test has the lowest number of healthy cows of the test 1 and 2 (Table 5.1 and Table 5.2). This test also has

(29)

Table 5.3: The distribution of healthy and ill cows in reality and the outcome of the algorithm for ﬁnding ill cows at test 3

The fourth and ﬁfth tests are done at the same farm (Tabels 5.4 and 5.5).

The fourth test has the highest specificity of all tests but a low sensitivity. This might however be occasional since this farm has the fewest cows. The fifth test has a higher sensitivity but a lower specificity.

Sensitivity 50%

Sensitivity 25%

The last test (Table 5.6) has the highest number of healthy cows and a high speciﬁcity. The sensitivity is however the lowest, but the number of ill cows are at this test also the lowest.

All these six tests results in a speciﬁcity of 77.0% and sensitivity 52.6%.

(30)

Chapter 6

Finding ill cows with a cusum-test

In this chapter it is presented how to detect cows with mastitis with a cusum- test. The basic idea is that if no teat has mastitis, the different teats will in average yield milk with the same conductivity, and they will in average yield with the same secretion rate. So each teat’s deviation from the average of the other teats conductivity and secretion rate will approach zero if it is summed up for each milking. If however one teat develops an inflammation, the conductivity will rise and the secretion rate will fall on that teat. Should this deviation from the average of the other teats be large enough, the inflammation can be detected.

6.1 Cusum test

Consider the test statistic c_t. How much the value diﬀers from v is summed up.

If c_tis equally distributed around v the sum will approach zero if not, s_t = c_t− v ,

g_t = g_t₋₁+ s_t ,

gt = 0 if gt< 0 , (6.1)

g_t = 0 and t_a= t and alarm if g_t> h > 0 .

So when the level h is reached at the time t_a, the change in the mean is detected.

Cusum test (6.1) could be used to detect changes in the ratios of the secretion rates and conductivity values.

6.2 The proposed algorithm

In the proposed algorithm the test statistic is the ratio of the diﬀerent teats conductivity (compared to each other), and the ratio of the diﬀerent teats’

(compared to each other) secretion rate. So if all teats have the same conductivity and secretion rate the ratio in all the four cusum-tests will be one. Hence, if one teat has mastitis the conductivity will ideally rise on that particular teat,

(31)

of how much the ratios diﬀer from one is a measurement on how ill the cow possibly is. When the deviation is large enough or occurs long enough an alarm is given, indicating the possibility of mastitis.

The proposed algorithm is based on a cusum-test. Parts from the current algorithm are used in the new algorithm as well. For example it is desired that a counter is included in the proposed algorithm. Also like in the current algorithm, the new method compares conductivity of the diﬀerent teats with each other. An idea was also to have a part of the ratio of secretion rate in the cusum test as well.

6.3 Modiﬁed cusum-test

For the detection of illness stis the ratio of the different teats secretion rate as well as conductivity, therefore it is necessary to have two inputs ctand ht. In (6.2), ct is the ratio between the different teats conductivity values, and ht is the ratio of the different teats secretion rate value. See Appendix A.1 for the calculations of the different ratios.

For each teat there will be a cusum test with, st = A(ct− vc) + B(ht− vh) , gt = gt−1+ st,

g_t = 0 if g_t< 0 , (6.2)

gt = 0 and ta= t and alarm if gt> h > 0 .

For the normal cusum-test the reset level h is chosen to 2. So when the sum has reached 2 the cow might be suffering from mastitis. In the modified cusum-test the sum is not immediately reset when the sum reaches 2, instead a counter is counted up. This leads to a modified cusum test with a counter.

6.3.1 Reset time

In a cusum-test the use of reset time is of great importance.

But in the detection of ill cows it is of less importance at what particular instant the cow is thought to be ill. Instead it is wished that an indication on the likelihood of illness is shown. The reset-time is normally the primary output in a cusum-test. In the proposed algorithm however the information noted is that the sum is over a certain level and the sum is not reset immediately. The critical level is really 2, but since the incorporation of a counter, the algorithm will continuously increase or decay until it reaches four, or zero. As soon as the sum is over two, the counter will be assigned a value between 0 and 10. If the sum is over four, the sum is reset. At what instant the sum is reset is however not of importance.

So s_t is summed up until the sum reaches a given level h (here 4), then the sum is reset to zero. A and B are positive constants. The drift terms vary slightly depending of if all teats are milked or not. If a conductivity value is below 3000µS and/or the secretionrate is below 50 g/h the teat is considered not milked, and ctand/or ht will be zero. If only one teat is milked there will be no comparison.

(32)

6.3.2 Parameters in the cusum-test

The parameters in the cusum-test were optimized to maximize the sum of sensitivity and specificity. Since in the current algorithm if one teat has a conductivity which is more than 15% higher than the average of the two lowest conductivity values, this was the first value for the parameter vc. But since too few cows were found vc was decreased. Because secretion rate vary more than conductivity, v_h will be larger. The parameters v_c and v_hwere optimized at two test so that the sum of the conductivity and secretion rate would be as large as possible, when doing this the parameters was A = 1 and B = 0 and then A = 0 and B = 1 for concluding v_h. Values from 1 increasing by 0.01 to 2 were used for the conclusion of the parameters v_h and v_c. So from the two tests the parameters v_c and v_h were decided. At three other test however the parameters A and B were optimized. Because of the difficulties in optimizing the parameters A, B, vcand vhat the same time other values of the parameters might give a slightly better result.

6.3.3 The calculations of conductivity ratio and secretion ratio

Depending on how many teats that are milked the ratios are calculated in diﬀer- ent ways. Also the drift terms vary slightly. If fewer than four teats are milked the ratios are a bit more unsure and the drift term needs to be larger as one has fewer data to compare between. The secretion-rate has a greater variance than conductivity. Therefore its drift term is larger. It is not unusual for a teat to have a secretion rate that is half the value of the others. The values for the drift terms and the ratios are found in Appendix A.1.

6.4 Results of the proposed algorithm

It was soon clear that it was not easy to reach the requirements of sensitivity of at least 80% and specificity of at least 90%. It is however shown that the results of the current algorithm is far from the requirements set by DeLaval as well. Tables 6.1-6.6 show that the specificity for the proposed algorithm varies between 78.3% and 93.3% and the sensitivity between 20.0% and 75.0%. The great variance between the smallest and highest sensitivity is due to test at farms where only few cows were ill (See Tables 6.4, 6.5 and 6.6). Tables 6.4 and 6.5 are tests from the same farm, but both have only a few ill cows. There is only one test done at the farm shown in Table 6.6. As stated in Section 5.2.3 results might vary depending of the conditions on the farm and breed kept at farm. By choosing A = 1 and B = 0, better results when it comes to specificity are gained and sensitivity is not much worse. If B > 0, both specificity and sensitivity will be lower than when B = 0. The fact that the secretion rate has a greater variance than the conductivity contributes to the difficulties to have secretion rate as a part of the cusum-test. The tests are the same as in Section 5.2.3.

Since the current and proposed algorithm are diﬀerent types of algorithms there will be cases where one algorithm is better than the other. Figures 6.1

(33)

There are also cases where the current algorithm is better than the proposed.

For example if an ill teat produce milk with occasionally high conductivity and occasionally low conductivity, the cow will be considered ill each time the teat produces milk with high conductivity (see Section 5.2.2 for explanation on how the current algorithm detects mastitis). If one teat has a conductivity that is for a long time slightly higher than the others’ the cusum-test might indicate that the cow is ill although it is not.

6.255 6.26 6.265 6.27 6.275 6.28 6.285 x 10⁴ 3600

3800 4000 4200 4400 4600 4800

Time [h]

RF LR RR

Figure 6.1: In this ﬁgure the conductivity is shown. This cow was ill on right- front teat (RF) when tested (the solid vertical line) and was found ill by the proposed but not the current algorithm.

These six test together results in a speciﬁcity of 85.4% and a sensitivity of 45.2%.

(34)

6.25 6.255 6.26 6.265 6.27 6.275 6.28 6.285 6.29 x 10⁴ 3500

4000 4500 5000 5500 6000

Time [h]

max(LF, RF and RR) LR

Figure 6.2: This graph shows conductivity. This cow was not ill when tested (the solid vertical line), which the proposed also found, the current algorithm however found the cow ill at teat left-rear (LR) The values of this teat is shown with the solid line. The stars represent the maximum values of the other teats conductivity.

Reality 92 Healthy Cows 13Ill Cows Outcome 72 Healthy 20 Ill 6 Ill 7 Healthy

(35)

Sensitivity 20%

Speciﬁcity 87%

Sensitivity 50%

1 2 3 4 5 6

0 20 40 60 80 100

Test nr

Procent

p−spec c−spec

1 2 3 4 5 6

0 20 40 60 80

Test nr

Procent

p−sens c−sens

Figure 6.3: Sensitivity and speciﬁcity of the current and proposed algorithm, c is current algorithm, p proposed, sens is sensitivty, spec is speciﬁcity

(36)

Chapter 7

Cows milked on only two or three teats

In this chapter it is investigated if it is possible to detect illness in cattle which are only milked on two or three teats. There is no counter in this test. This algorithm takes both conductivity and secretion rate into consideration. It was previously thought that the current algorithm was not as accurate for cows milked on less than four teats as it was for cows milked on four teats. It is in this chapter however shown that the current algorithm has a high sensitivity and speciﬁcity on cows milked on less than four teats. Because it was thought that the current algorithm had a poor outcome it should be investigated if a special algorithm for cows milked on less than four teats could be developed.

7.1 Proposed algorithm

In this algorithm secretion rate and conductivity are weighted together to ﬁnd mastitis. For all milkings in a time interval points are set for secretion rate and conductivity, these points are added up and divided by the number of milkings in the interval. If the value then exceeds a certain value the cow is considered ill.

The points represent how much conductivity and secretion rate are weighted. If a cow is milked infrequently and has low secretion rate and high conductivity at one milking it has a greater impact on the algorithm. The points are diﬀerent depending on if it is secretion rate or conductivity. Conductivity is weighted more when it comes to ﬁnding mastitis. So the points represent how much high conductivity and low secretion rate are weighted.

7.1.1 Time interval for measurements

When developing this algorithm the interval in which the values of the conductivity and the secretion rate are used for determining if the cow is ill or not, is

±24hours of the test date.

(37)

7.1.2 Points

The points (or weights) are optimized so that the sum of sensitivity and speci- ﬁcity will be as large as possible. The points are optimized at three tests and veriﬁed at nine tests.

If only two teats are milked and the largest value of conductivity is 1.15 times the smallest, one point is set for the teat with highest conductivity. If three teats are milked and the highest conductivity is 10% higher than the smallest one, one point is set. If the conductivity excedes 7000µS/cm for any of the teats two points are set. If the teat with the lowest secretion rate is 30%

lower than the highest secretion rate, another 0.5 points is added to the points received above. The teat with 30% higher secretion rate than the second highest secretion rate gets 0.5 points subtracted if three teats are milked. If only two teats are milked the value has to be 70% higher than the second highest to get 0.5 points subtracted. If the sum after it has been divided with the number of milkings done in 48h is over 0.9 the cow is considered possibly ill.

7.2 Results

Table 7.1 shows how many cows which in the interval are milked on only two or three teats (# cows) and how many of those which are ill (# ill cows). Healthy cows pointed out as healthy (healthy) are also shown, as how many of the ill cows pointed out as ill (ill found). The same applies in Table 7.2 which shows the outcome of the current algorithm at these tests. The current algorithm for cows milked on less than four teats is the same as described in section 5.2.2.

Farm Date # cows # ill cows ill found healthy

Farm 1 Test 1 20 4 4 0

Farm 1 Test 2 22 3 2 2

Farm 1 Test 3 13 3 2 0

Farm 1 Test 4 18 4 3 2

Farm 1 Test 5 5 2 1 2

Farm 2 Test 6 9 3 3 0

Farm 2 Test 7 7 1 0 0

Farm 3 Test 8 20 2 1 1

Farm 4 Test 9 5 0 0 2

total 119 22 16 9

Table 7.1: The outcome of the proposed algorithm

That means that the algorithm found 72.7% of the ill cows. The algorithm points out 9.3% of the healthy cows as ill.

So from Table 7.1 and Table 7.2 it is clear that both algorithms have a sensitivity of 72.7% whereas the proposed algorithm has a speciﬁcity of 90.7%

and the current has a speciﬁcity of 89.7%.

(38)

Farm Date nof cows nof ill cows ill found healthy

Farm 1 test 1 20 4 3 3

Farm 1 test 2 22 3 2 3

Farm 1 test 3 13 3 1 0

Farm 1 test 4 18 4 4 1

Farm 1 test 5 5 2 1 2

Farm 2 test 6 9 3 3 0

Farm 2 test 7 7 1 1 0

Farm 3 test 8 20 2 1 1

Farm 4 test 9 5 0 0 0

total 119 22 16 10

Table 7.2: The outcome of the current algorithm

(39)

Chapter 8

Conclusions on ﬁnding ill cows

Mastitis is an illness which is characterized through the inflammation of the mammary gland. The current algorithm considers a teat as a possible mastitis carrier if its milk has a high conductivity value or has a relative high conductivity value when compared to the other teats’. The current algorithm has a lower sensitivity and specificity than preferable. This thesis work investigates if it is possible to also take secretion rate into consideration when concluding if a teat might have mastitis. Since it is thought that the current algorithm has a particular low sensitivity and specificity on cows milked on fewer than four teats, it is also investigated if a special algorithm for cows milked on fewer than four teats could be developed.

A cusum test could be used to detect slowly changes in the ratios of the different teats conductivity values and secretion rate. It is, although the current algorithm for finding ill cows is not very accurate, hard to find a better algorithm. An algorithm based on a cusum test with only a few percent higher specificity has been presented. The current algorithm has an average (at the six test conducted) of specificity 78.2% and sensitivity 45.2% whereas the proposed algorithm has specificity 86.1% and sensitivity 47.3%. The proposed algorithm has a slightly larger sensitivity because of one test with few cows where the proposed algorithm found more of the ill cows. When looking at all the cows together the current algorithm has a specificity of 77.0% and a sensitivity of 51.7%. The proposed algorithm gets a specificity of 85.4% and a sensitivity of 45.2%. So when summing up specificity and sensitivity, the proposed algorithm is only slightly better. This shows the difficulties in developing a new algorithm for the detection of ill cows. The algorithm specially developed for cows milked on fewer than four teats, also shows only slightly improvements when compared to the existing algorithm.

One could also draw the conclusion from the algorithm for cows milked on less than four teats, that taking secretion rate into consideration does not improve algorithms for ﬁnding mastitis a great deal.

The modiﬁed cusum algorithm gets an outcome that varies a great deal at different farms. It is however clear from the proposed algorithm for expected yield

(40)

that diﬀerent farms get diﬀerent results. There is it shown that in how many cases the proposed algorithm is better varies with as much as 24.6%. If more tests were conducted a better understanding on why the algorithm is better at some farms would be gained. Also if the constants in the cusum-test were optimized from more tests a better overall result might be the outcome.

(41)

Bibliography

[1] Eﬃcient Milking, Kerstin Svennersten-Sjaunja, DeLaval 2001.

[2] Discrete-time Stochastic Systems Estimation and Control, Torsten S¨oder- str¨om, Prentice Hall 1994.

[3] Machine Milking and Lactation, editor A. John Bramley et al, Jorn Ham- man & Frank H. Dodd, Insight Books 1992.

[4] Adaptive Filtering and Change Detection, Fredrik Gustafsson , John Wiley

& Sons LTD 2001.

(42)

Appendix A

Appendix

A.1 Drift terms and ratios for the modiﬁed cusum- test

A.1.1 Three teats milked

If three teats are milked, the drift vc = 1.12 and vh= 1.35 .

If teat Z has the highest or second highest conductivity value c^Z_t = 2∗ conductance value of teat Z smallest value+second smallest. If teat Z has the lowest conductivity value c^Z_t =conductance value of teat Z

highest value

If Z has the lowest secretion rate h^Z_t =highest secretion rate+second highest 2∗secretion rate teat Z . If Z has the highest secretion rate h^Z_t = lowest secretion rate+second lowest

2∗secretion rate teat Z

A.1.2 Two teats milked

If two teats are milked, the drift vc = 1.12 and vh= 1.35 .

If teat Y has the highest conductivity value c^Y_t = 2∗conductance value of teatY smallest value . If Y has the lowest conductivity value c^Y_t =conductance value of teatY

highest value . If Y has the lowest secretion rate h^Y_t = highest secretion rate

secretion rate of teat X. If Y has the highest secretion rate h^Y_t = lowest secretion rate

secretion rate of teat X.

(43)

A.1.3 Four teats milked

If four teats are milked, the drifts will be vc= 1.1 and vh= 1.3 The ratios will be as follows:

If teat X has the highest or second highest conductivity value c^X_t = 2∗

conductance value of teat X smallest value+second smalles value.

If however teat X has the lowest or second lowest conductivity value. c^X_t = 2∗ conductance value of teatX

highest value+second highest value

If teat X has the lowest or second lowest secretion rate h^X_t = highest secretion rate+second highest 2∗secretion rate of teat X . If teat X has the highest or second highest secretion rate h^X_t =lowest secretion rate+second lowest

2∗secretion rate of teat X .

On Predicting Milk Yield and Detection of Ill Cows