The SimBatch simulation model Anders Ringgaard Kristensen - - PowerPoint PPT Presentation

The SimBatch simulation model

Anders Ringgaard Kristensen

Department of Large Animal Sciences

System description

• A batch of slaughter pigs is simulated

from insertion of the batch until slaughter.

• The pigs are fed ad libitum.
• Sent to a Danish slaughter house.
• The highest price per kg is obtained in

the “optimal” interval from 70 to 84 kg slaughter weight.

• Delivery is decided on the basis of
• bserved live weight.
• The slaughter policy is decided by the

farmer

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 2

Model type SimBatch is a

• Dynamic (i.e. not static) time stepping model

with daily updating of the states of the pigs and the batch.

• Mechanistic (i.e. not empirical) model since it

models a batch by modeling each individual pig.

• Initially it is a deterministic model, but in the

exercise it is step by step transformed into a stochastic model.

• Monte Carlo simulation model, because the

simulation mechanism is based on drawing random numbers from relevant distributions.

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 3

State of nature, Φ0

# State of nature information stateOfNature = list(initialAverage = 30, # Average weight at insertion initialStdDev = 2, # Std. deviation of weight at insertion averageDG = 0.900, # Average daily gain stdDevDG = 0.100, # Std. dev. of daily gain between pigs autocorrelationDG = 0.99, # Autocor. for individual daily gain diseaseRisk = 0.01, # Risk of disease any day diseaseEffectAve = 0.250, # Average effect of disease diseaseEffectStd = 0.030, # Standard deviation of disease effect diseaseDurationAv = 5, # Average duration of disease liveToSlaughterAv = 0.76, # Live w. to slaughter w. conversion liveToSlaughterStd = 1.4, # Standard deviation of conversion lowerOptimal = 70, # Lower bound of opt. slaughter weight upperOptimal = 84, # Upper bound of opt. slaughter weight deliverEvery = 7) # Days btw. del.(7 = weekly, 1 = daily)

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 4

Decision strategy, Θ

# Decision strategy information decisionStrategy = list(thresholdWeight = 100, # Intended live weight at slaughter minimumDelivery = 10, # Minimum number of pigs to send to maximumUnderWeight = 2, # Acc this number of pigs bel. Thresh. minimumStock = 50, # Minimum acc. stock: Send all if below maximumAge = 90) # Maximum age before slaughter # Strategy with completely individual delivery individualStrategy = list(thresholdWeight = 100, # Intended live weight at slaughter minimumDelivery = 1, # Minimum number of pigs to send to maximumUnderWeight = 0, # Acc this number of pigs bel. thres minimumStock = 0, # Minimum acc. stock: Send all if bel maximumAge = 9999) # Maximum age before slaughter # Strategy with slaughter at a given age ageStrategy = list(thresholdWeight = 9999, # Intended live weight at slaughter minimumDelivery = 9999, # Minimum number of pigs to send to maximumUnderWeight = 0, # Acc this number of pigs bel. thres minimumStock = 0, # Minimum acc. stock: Send all if bel maximumAge = 90) # Maximum age before slaughter

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 5

Modeling a pig

A pig is defined by 4 state variables

• Live weight in kg
• The gain today in kg (if healthy)
• An integer, i, defining the status of the pig:
• i = -2: The pig has already been slaughtered
• i = -1: The pig is healthy and still present
• i ≥ 0: The pig is diseased and will stay diseased i

days yet

• The effect of disease (if any) on daily gain

In SimBatch, a pig is represented on a given day by a vector with 4 elements. The state of the pig is updated every day.

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 6

Creating and printing a pig A valid pig object: > pig = c(45.0, 0.956, 3, 0.123) > pig [1] 45.000 0.956 3.000 0.123 A pig with a live weight of 45 kg, a gain until tomorrow of 956 g (if healthy). It is, however diseased and will stay diseased for 3 days. The effect

• f disease is a reduction in daily gain of 123 g

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 7

Simple functions used with a pig object

getStartWeight(stateOfNature): A start weight (in kg) is returned. Initially, it only returns the initial average from the state of nature. getInitialDG(stateOfNature): An initial daily gain (in kg) for a pig is

• returned. Initially, it only returns the average daily gain from the state of

nature. isPigDiseased(stateOfNature): Returns a logical value. If the pig gets diseased it returns TRUE, otherwise FALSE. Initially, it always returns FALSE. getDiseaseEffect(stateOfNature): Returns the reduction (in kg) of daily gain during a disease period. Initially it only returns the average disease effect from the state of nature. getDiseaseDuration(stateOfNature): Returns the duration (in days) of a new disease. Initially, it only returns the average duration from the state of nature. getSlaughterWeight(liveWeight, stateOfNature): Transforms the specified input parameter liveWeight to slaughter weight (in kg). Initially it just multiplies the live weight with the conversion ratio from the state of nature. getObservedWeight(liveWeight): Returns the observed live weight of a pig with the specified true live weight. Initially, it just returns the true live weight (ignoring the observation error). updateDG(oldDG, stateOfNature): Returns a new daily gain for the next 24 hours from the specified “old” daily gain (from the previous 24 hours). Initially, it just returns the old value (without any transform).

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 8

Function creating a piglet (don’t change)

# Create a "30 kg" piglet at random createPiglet = function(son) { piglet = rep(NA, 4) piglet[1] = getStartWeight(son) # Draw initial live weight piglet[2] = getInitialDG(son) # Draw initial daily gain piglet[3] = -1 # Number of days remaining if diseased # If "-1": Healthy. If "-2": Slaughtered piglet[4] = 0 # Disease effect on daily gain return(piglet) }

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 9

Function updating a pig (don’t change)

# Updates the state of a pig updatePig = function(pig, son) { pigNew = rep(NA, 4) If (pig[3] > -2) { # Is the pig still present? pigNew[1] = pig[1] + pig[2] - pig[4] # Yes, calculate new weight pigNew[2] = updateDG(pig[2], son) # Update daily gain if (pig[3] == -1) { # Is the pig healthy now? if (isPigDiseased(son)) { # Yes it was. Should it change state? pigNew[3] = getDiseaseDuration(son) # Yes it should. Draw duration of disease pigNew[4] = getDiseaseEffect(son) # Draw effect of disease } else { # No, the pig remains healthy pigNew[3] = -1 # Define it as healthy pigNew[4] = 0 # No effect of disease } } else { # No, the pig is already diseased if (pig[3] > 0) { # Has it still at least one day left of disease? pigNew[3] = pig[3] - 1; # Yes, reduce the number of days left by one pigNew[4] = pig[4] # Keep the disease effect on daily gain } else { # No, it is the last day of the disease pigNew[3] = -1; # Change state to healthy pigNew[4] = 0 # No disease effect next time } } } else { # The pig is already slaughtered pigNew = pig # Leave it as it is } return(pigNew) }

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 10

Representation of a batch on a given day

Just an array of pigs. Example of a batch with 5 piglets:

> pig1 = createPiglet(stateOfNature) > pig2 = createPiglet(stateOfNature) > pig3 = createPiglet(stateOfNature) > pig4 = createPiglet(stateOfNature) > pig5 = createPiglet(stateOfNature) > batchDay = array(NA, dim =c(5, length(pig1))) > batchDay[1,] = pig1 > batchDay[2,] = pig2 > batchDay[3,] = pig3 > batchDay[4,] = pig4 > batchDay[5,] = pig5 > batchDay [,1] [,2] [,3] [,4] [1,] 29.77745 0.8219607 -1 0 [2,] 27.78206 0.9373041 -1 0 [3,] 27.52593 0.8217774 -1 0 [4,] 31.96445 0.9472012 -1 0 [5,] 30.67659 0.8425552 -1 0

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 11

Function updating a batch (don’t change)

# Updates the state of a batch updateBatch = function(batchDay, son) { newBatchDay = array(NA, dim = c(length(batchDay[,1]), length(batchDay[1,]))) for (i in 1:length(batchDay[,1])) { # Iterate over pigs pigDay = batchDay[i,] # Take the i’th pig pigDay = updatePig(pigDay, son) # Update it newBatchDay[i,] = pigDay # Insert it in the updated batch } return(newBatchDay) # Return the updated batch }

Function counting pigs left in a batch (don’t change)

# Counts the number of pigs still present in the batch (not slaughtered) pigsLeft = function(batchDay) { left = 0 for (i in 1:length(batchDay[,1])) { if (batchDay[i,3] > -2) { left = left + 1 } } return(left) } Department of Large Animal Sciences

Slide 12

Applying decision strategy to batch on a given day

# Applies the decision strategy on the batch a given day selectForSlaughter = function(batchDay, age, dec) { correctedDay = array(NA, dim=c(length(batchDay[,1]), length(batchDay[1,])))

• bservedDay = rep(NA, length(batchDay[,1]))

for (i in 1:length(batchDay[,1])) { correctedDay[i,] = batchDay[i,] # Copy pigs to new array

• bservedDay[i] = getObservedWeight(batchDay[i,1])

# Get observed weights } leftPigs = pigsLeft(batchDay) # How many pigs are still present if (leftPigs < dec$minimumStock || # Fewer than minimum? age > dec$maximumAge) { # Maximum age exceeded? for (i in 1:length(batchDay[,1])) { # Yes, correctedDay[i,3] = -2 # ... cull the rest } } readyIndexes = which(observedDay >= dec$thresholdWeight) # Find the pigs with observed weights exceeding threshold readyCount = length(readyIndexes) # How many were there if (readyCount < dec$minimumDelivery) { # Are there two few? if (readyCount > dec$minimumDelivery - dec$maximumUnderWeight) { # Yes, but are there so many that we slaughter anyway? sortedVector = sort(observedDay, decreasing = TRUE) # Yes, find the biggest pigs realThres = sortedVector[dec$minimumDelivery] # Find the real threshold readyIndexes = which(observedDay >= realThres) # Pick all pigs exceeding the real threshold for (i in readyIndexes) { correctedDay[i,3] = -2 # Mark them as slaughtered } } } else { # No, there were enough pigs for (i in readyIndexes) { # Pick them correctedDay[i,3] = -2 # Mark them as slaughtered } } return(correctedDay) # Return the pigs with the slaughtered marked }

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 13

A simulation run – function simulateBatch (don’t change)

# Simulates a batch with the specified number of pigs over the specified number of days simulateBatch = function(pigs, days, son, dec) { piglet = createPiglet(son) # Create a piglet (just to know its dimension ...) batchStructure = array(NA, dim=c(pigs, days, length(piglet))) # Create a data structure for the results for (i in 1:pigs) { # Create "pigs" piglets piglet = createPiglet(son) # The i’th piglet batchStructure[i,1,] = piglet # Insert it } for (n in 2:days) { # Iterate over the simulation period day by day batchStructure[,n,] = updateBatch(batchStructure[,n-1,], son) # if (n %% son$deliverEvery == 0) { # Is it a delivery day? batchStructure[,n,] = selectForSlaughter(batchStructure[,n,], n, dec) # Yes, apply decision strategy } } return(batchStructure) # Return full data structure }

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 14

Plotting weight curves from a simulation run

# Plot growth curves of pigs plotBatchWeights = function(batch, skip) { pigWeights = batch[1,,] days = 1:length(pigWeights[,1]) plot(days, pigWeights[,1], type="l", ylim=c(25, 110), col="white", xlab="Days", ylab="Kg live weight") for (i in 1:length(batch[,1,1])) { pigLife = batch[i,,] liveIndexes = which(pigLife[,3] > -2) pigWeights = pigLife[liveIndexes,1] if ((i-1) %% skip == 0) { lines(pigWeights) } } }

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 15

Plottet weight curves from a simulation run

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 16

Real data (automatic weighing system) - Stygar et al. PigIT

Department of Large Animal Sciences

Advanced Quantitative Methods in Herd Management Slide 17

Time in hours

Functions returning output variables

batch = simulateBatch(100, 120, stateOfNature, decisionStrategy)

• Runs a simulation with 100 pigs over 120 days using the specified

state of nature and the specified decision strategy. plotBatchWeights(batch, 1)

• Plots the resulting weight curves

deliveredInOptimalInterval(batch, stateOfNature)

• Calculate (and print) the number of pigs delivered in the optimal

interval dailyGain(batch)

• Calculates the true average daily gain for the batch

Advanced Quantitative Methods in Herd Management

Department of Large Animal Sciences

Slide 18