CS/ECE/ISyE 524 — Introduction to Optimization — Spring 2016

Minimizing Economic and Public Health Impacts of Radioactive Contamination

AJ Nosek ([email protected]) and Newton Wolfe ([email protected])


Table of Contents

  1. Introduction
  2. Mathematical Model
    1. Cost of Decontamination Methods
    2. Dose Contributions from Contaminated Surfaces
    3. Representative Areas Based on Population
    4. Model
  3. Solution
    1. Decontaminating Fukushima
    2. Pareto Tradeoff
  4. Results and Discussion
  5. Conclusion

1. Introduction

While an accident at nuclear power plant leading to an uncontrolled release of radioactive material into the atmosphere is an unlikely event, it can also cause extensive consequences.  Our project is to create an optimization tool that allows a user to minimize offsite cleanup costs caused from a nuclear power plant accident, such as the Fukushima nuclear disaster.  A number of different decontamination methods exist for various types of surfaces and with different characteristics such as effectiveness, costs, and constraints.  Our tool allows a user to enter the amount of contaminated area that exists for their situation (broken down by land usage), and the amount of cleanup dose reduction necessary to reach acceptable dose levels.  Our tool will then recommend decontamination methods for the different areas that will minimize cleanup costs.  An example is provided using GIS information of the contaminated areas around Fukushima detailing the contamination levels and different land use areas.  

Most planning for radiological accidents focus on the emergency response, that is, making sure the public is adequately protected from radiation.  However, arguably the more consequential aspect of a nuclear power plant accident is the long-term issue of offsite contamination.  Contamination can displace individuals from their homes and livelihoods, and this disruption can be traumatic and expensive.  

About 150,000 people evacuated from Fukushima, including about 76,000 people from areas that were not allowed to return because their homes are in areas that have been restricted.  Five years after the accident has occurred, only a fraction of these areas have had their restrictions lifted.  Arrangement of the restricted areas and decontamination plans of the various contaminated provinces took months to years to finish, and decontamination is still ongoing in many places.  All the while, many displaced individuals have been living in limbo.  Being displaced and living in temporary housing has separated families and livelihoods. [Japan Ministry of the Environment]

These conditions show the need for better recovery planning and response.  By creating a customizable tool that can minimize decontamination costs, decision makers can make better informed decisions for their situation regarding how to best decontaminate areas.  It will also help them decide when to condemn areas that may not economical or feasible to restore when necessary, and instead help people restore their lives elsewhere.  

2. Mathematical model

2.A. Cost of Decontamination Methods

In order to estimate the decontamination costs, first we had to identify different decontamination methods and estimate their costs. The list of decontamination methods we used can be found in Andersson, et al (2003). The decontamination methods are specific to a type of surface. For instance, deep plowing a field is a method for open areas, replacing a roof is specific to roofs, and so forth. Each method has its own cost, effectiveness, and time period of applicability:

  • The time period of applicability is a constraint. For instance, cutting the grass is only effective for about the first week. The more time that passes, the radioactive aerosols responsible for long-term contamination (i.e. Cs-134/137) permeate the soil and other methods will need to be considered. Other methods may have similar time restrictions
  • The effectiveness of a decontamination method is known as the decontamination factor (DF). The decontamination factor will reduce the dose contributions from that surface, but this needs to be weighted by the importance of this surface compared to all the other contaminated surfaces in the area.
  • The cost of a decontamination method is the sum of the equipment, labor, and fuel costs. The labor costs are based on the average salary for a Hazardous Materials Removal Worker, as reported by the Bureau of Labor Statistics. Labor wages are doubled to account for overhead. It is assumed that workers work 10 hour days, 5 days a week working, and spend 60% of their time performing remediation activities, while the other 40% is spent in travel into the remediation zone, donning protective gear, and miscellaneous activities. The equipment costs are based on the type of decontamination (e.g. road planer), and are normalized into a daily cost (assuming a decontamination period of no fewer than 0.3 years). Fuel costs are based on the type of equipment being used. All costs are considered as a function of the surface decontamination area (\$/$\textrm{m}^2$), by first computing the hourly cost (\$/hr) and then multiplying by the work rate of the specific decontamination method (h/$\textrm{m}^2$). All costs are informed by Anderrson (2003), internet queries, and personal judgment. The user can update these costs are appropriate in the attached spreadsheet: “NosekWolfe_countermeasures.xls”

2.B. Dose Contributions from Contaminated Surfaces

The dose received by an individual will depend on two things: the distribution of contamination in the environment, and the importance of different surfaces. The distribution of contamination is a function of the relative deposition from the plume of radioactive material. Trees and soils will be relatively more contaminated than inside a housing unit. However, the interior of a housing unit will have closer proximity to individuals. The contributions from different surfaces were calculated using MCNP, a particle transport code, using a number models to represent different housing environments, such as wood house, a brick house, a row of terrace housing units, and an apartment complex. The dose contributions from different surfaces are based on Andersson (1996), and were updated to include doses from interior surfaces based on Nisbet et al (2007). We also estimated a set of dose contributions based on these other models and judgment for a farm house property. The relative dose contributions from contaminated surfaces in the representative properties are below:

Contaminated Surfaces: Roads Residential Soil Open Soil Trees Walls Roofs Interior
Farm House 0.00 0.15 0.45 0.05 0.08 0.11 0.16
Wood House 0.00 0.39 0.00 0.29 0.07 0.09 0.16
Brick House 0.00 0.39 0.00 0.28 0.04 0.13 0.16
Semi-detached House 0.00 0.45 0.00 0.22 0.03 0.08 0.16
Terrace (town) House 0.09 0.42 0.00 0.22 0.02 0.16 0.16
Apartment Building 0.18 0.42 0.00 0.20 0.03 0.00 0.16

2.C. Representative Areas Based on Population

A challenge in formulating this program was relating contaminated geographic areas to local contaminated surfaces that contribute to dose. A number of different geographic areas were modeled based on population density. These geographic areas define a characteristic model, largely based on housing density. Seven different characteristic areas were developed for increasing population densities. For instance, the lowest density area considers a wood house with a standard footprint on 10 acres. From this characteristic model, the areas of all surfaces are estimated. This property is used to fill the 50-100 people / $\textrm{km}^2$ density. More dense areas use smaller property areas, and higher density housing units.

The following table shows the estimated amount of surface area in a one square kilometer area for the different characteristic population densities:

Area Type Population Density (/km^2) Roads (m^2) Residential Soil Open Soil Trees Walls Roofs Interior
farm-1 50-100 0 75000 925000 0 0 0 0
low-2 100 – 200 76168 456184 456184 0 10329 11465 52779
low-3 200 – 500 169261 805262 0 402631 22953 25477 117287
low-4 500 – 1000 179218 766830 0 383415 48605 53952 248373
med-5 1000 – 2000 158723 733373 0 366686 97211 107904 496747
med-6 2000 – 3000 269830 546734 0 273367 165258 183437 844469
high-7 3000+ 197677 153684 0 76842 495775 648638 2631737

2.D. Model

The model is an integer linear program because only one decontamination method may be used on a surface in a given area, as the entire area must have the dose contribution reduced evenly to ensure that everyone's does is appropriately reduced.

The objective is to minimize overall costs incurred via the different decontamination methods. Since these costs are measured in a per square meter basis, we must scale it by the number of square meters of the surface being decontaminated per square kilometer of land to be decontaminated.

Similarly, the decontamination factors for each method must be scaled by the relative contribution of the surface that is being decontaminated to a person's overall dose, and the sum of all of these scaled decontamination factors must be at least as large as the required overall dose reduction.

$$ \begin{aligned} \underset{x}{\text{minimize}}\qquad& \sum p_irx_i && \; \; \; \qquad \textrm{where $p_i$ is the cost of decontamination method $x_i$ per square meter} \\ & && \; \; \; \qquad \textrm{and $r$ is the area in square meters that the surface covers per square kilometer } \\ \text{subject to:}\qquad& x_i \leq a_i && \forall i \qquad \textrm{where $a_i = 1$ if $x_i$ is allowed}\\ & gx_i \leq t_i && \forall i \qquad \textrm{where $t_i$ is the latest time that $x_i$ can be used} \\ & && \; \; \; \qquad \textrm{and $g$ is the time that decontamination can begin} \\ & \sum_{i\in S} x_i \leq 1 && \; \; \; \qquad \textrm{for every surface $S$} \\ & \sum cd_ix_i \geq t && \; \; \; \qquad \textrm{where $c$ is the surfaces's contribution to dose,} \\ & && \; \; \; \qquad \textrm{$d_i$ is the surface decontamination factor of $x_i$,} \\ & && \; \; \; \qquad \textrm{and $t$ is the required dose reduction} \\ & x_i \in \{ 0,1 \} && \forall i \end{aligned} $$

3. Solution

Because every decontamination scenario is different, our model allows the user to customize all of the setup. In this instance, we've extracted data from the Fukushima incident, based on radiation distribution from the Nuclear Regulation Authority and extrapolated land use data from the 2005 Census population map converted via analogy to researched dose exposures based on population density.

This part of the code also sets up the surfaces and decontamination methods we'll be using in this instance. The surfaces are somewhat fixed, in that the dose contribution and surface area models rely upon them, but the optimization itself is mostly setup-agnostic.

In [1]:
# Data about dose contributions and areas from models that we pull in for the optimization.

# Dose Contributions from contaminated surfaces (dry deposition)
# Anderrson (1996), updated with interior doses using data from EURANOS 2007.
# roads, rsoil, osoil, trees, walls, roofs, interior
doseUN = [ 0.00 1.00 0.00 0.00 0.00 0.00 0.00 ] # Uninhabited areas
doseFM = [ 0.00 0.15 0.45 0.05 0.08 0.11 0.16 ] # Farm (open field)
doseWH = [ 0.00 0.39 0.00 0.29 0.07 0.09 0.16 ] # Wood house
doseBH = [ 0.00 0.39 0.00 0.28 0.04 0.13 0.16 ] # Brick house 
doseSD = [ 0.00 0.45 0.00 0.22 0.03 0.08 0.16 ] # Semi-detached house
doseTH = [ 0.09 0.42 0.00 0.22 0.02 0.16 0.16 ] # Terrace (Town) house
doseAP = [ 0.18 0.42 0.00 0.20 0.03 0.00 0.16 ] # Apartment building

# Surface areas of characteristic model locations (NUREG/CR-xxxx), 
# expanded for Fukushima population densities
# roads, rsoil, osoil, trees, walls, roofs, interior
AreaUN0 = [0 0 0 0 0 0 0] # Uninhabited areas
AreaFM1 = [0 75000 925000 0 0 0 0] # Farm (50-100 people / km2)
AreaWH2 = [76168 456184 456184 0 10329 11465 52779] # Wood House (100-200 people / km2)
AreaWH3 = [169261 805262 0 402631 22953 25477 117287] # Wood House (200-500 people / km2)
AreaWH4 = [179218 766830 0 383415 48605 53952 248373] # Wood House (500-1000 people / km2)
AreaTH5 = [158723 733373 0 366686 97211 107904 496747] # Terrace (Town) house (1000-2000 people / km2)
AreaTH6 = [269830 546734 0 273367 165258 183437 844469] # Terrace (Town) house (2000-3000 people / km2)
AreaAP7 = [197677 153684 0 76842 495775 648638 2631737] # Apartment Building (3000+ People / km2)

# Read in input file containing the decontamination methods
deconFile = readdlm("NosekWolfe_countermeasures.csv",',');
(m,n) = size(deconFile)

# Creates an array of distinct surfaces
# This is actually fixed and the above defined arrays depend upon it, but this ensures compatibility
# in our later searches.
surfaces = []
for i in 2:m
    if (findfirst(surfaces,deconFile[i,2]) == 0)
        push!(surfaces,deconFile[i,2])        
    end
end

n_method = 2:m      # rows containing methods
n_roads = 2:4       # rows containing methods for roads
n_rsoil = 5:13      # rows containing methods for residential soil/grass
n_osoil = 14:17     # rows containing methods for open soil/grass areas
n_trees = 18:19     # rows containing methods for trees
n_walls = 20:23     # rows containing methods for exterior walls
n_roofs = 24:27     # rows containing methods for roofs
n_inside = 28       # rows containing methods for interior surfaces

methodName = deconFile[n_method,1] # Name of the decontamination method
methodSurf = deconFile[n_method,2] # Surface to which the decontamination method applies
methodDF   = convert(Array{Float64,1},deconFile[n_method,3]) # Cast decontamination factors to floats
methodCost = convert(Array{Float64,1},deconFile[n_method,4]) # Cast costs to floats
methodTime = convert(Array{Float64,1},deconFile[n_method,5]) # Cast times to floats
# We'll also need an array corresponding to which methods the user allows.
methodUse  = ones(length(n_method),1);

The following blocks are variables that are easy for novice users to set because they don't require in-depth knowledge of decontamination procedures or radiation modeling. You can see how even just small changes can make big differences in costs, and in some cases can make it infeasible for humans decontaminate an area to a high degree. (Luckily, natural processes like weathering and radioisotope decay also reduce the dose as well. Furthermore, there are far more decontamination possibilities than could fit into a model like this; the few dozen that are imported above are a small subset of the possiblities.)

In [2]:
# The following variables are user-configurable even by novice users:

# Setting methodUse[i]=0 to disallow the ith countermeasure.
# For example, this disallows snow removal
methodUse[9] = 0

# The number of months after the radiation release
# Dictates which decontaminations measures can and can't be used
givenTime = 0.5;

The following block is the optimization routine. It takes in the amount of various surfaces in 1 square kilometer and the amount that each surface contributes to the average person's dose living in that area, and finds the minimum cost to decontaminate to reach a minimum dose reduction. It will also throw infeasibility if the required dose reduction is impossible.

The value it outputs is the minimum cost to decontaminate 1 square kilometer of land with the surfaces and dose contributions given to the degree required by the dose reduction given.

In [3]:
using JuMP
function solveOpt(area, doseContrib, doseReduction)
    m = Model()

    # Set 1 binary variable per each method in each surface
    @variable(m, use[1:length(n_method)], Bin)

    # Use at most one method per surface. (We can use zero if other surfaces make up for it.)
    @constraint(m, sum{use[i], i in (n_roads)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_rsoil)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_osoil)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_trees)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_walls)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_roofs)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_inside)-1} <= 1)

    # Limit methods to ones the user allows
    @constraint(m, use .<= methodUse) 
                            
    # Limit methods based on time.
    for i in 1:length(n_method)
        @constraint(m, use[i]*givenTime <= methodTime[i])
    end
    
    # Must decontaminate by at least the required decontamination factor
    @constraint(m, sum{use[i]*methodDF[i]*doseContrib[findfirst(surfaces,methodSurf[i])], i in 1:length(n_method)} >= doseReduction)

    # Want to minimize the cost based on cost of decontaminating every surface
    @objective(m, Min, sum{use[i]*methodCost[i]*area[findfirst(surfaces,methodSurf[i])], i in 1:length(n_method)})
    status = solve(m)
    return getobjectivevalue(m)
end;

3.A. Decontaminating Fukushima

In this block, we pulled in geographic data (through GIS software) of the Fukushima radiation release, and compared the population densities with where radiation was released. This data is then fed to the models seen in an earlier code block that determine how much of the dose each surface contributes, which is used by the optimization routine to find the minimum-cost subset of decontamination methods.

This is a map of the radiation data we used: Radiation distribution map And this is a map of population density data we used overlayed on top: Radiation distribution map with overlayed population data

In the end, this block then returns the total cost to decontaminate the entire area.

In [4]:
# Now, we get to bring it all together!
# This input matrix is user-customizable. Here we use data from Fukushima circa April 2011
# Note that it's easy to set an overly-agressive decontamination factor
# If you do this you'll get an error from the solver and no total cost
# Unfortunately, JuMP and Julia don't appear to have a way to handle them

# The first row of this matrix is the required dose reduction. Each row after that is a type of populated
# area whose dose may need reducing, and each cell represents a certain area with a certain dose.
# The contents of the cell are the number of square kilometers of that dose.
input = [
0.8         0.8         0.6         0           0           0           0           0
0           0           0           0           167.586696  100.5520176 33.5173392  33.5173392
2.62275819  2.62275819  7.86827457  7.86827457  104.9103276 31.47309828 52.4551638  52.4551638
4.25753427  8.51506854  8.51506854  21.28767135 85.1506854  127.7260281 85.1506854  85.1506854
12.55895075 12.55895075 37.67685225 62.79475375 142.34991   350.2501371 323.7317517 313.9737688
2.62275819  5.24551638  7.86827457  18.35930733 39.34137285 44.58688923 78.6827457  65.56895475
12.67847151 12.67847151 12.67847151 38.03541453 177.4986011 253.5694302 507.1388604 253.5694302
140.4190676 84.25144056 196.5866946 280.8381352 140.4190676 421.2572028 1123.352541 421.25720282
]

# These two map population densities to the residential areas which they represent
# Residential areas are the vast majority of the development near Fukushima
areamapping = Array[ AreaTH6, AreaTH5, AreaWH4, AreaWH3, AreaWH2, AreaFM1, AreaUN0]
# This dose mapping maps the areas above to the appropriate dose contributions
dosemapping = Array[ doseTH,  doseTH,  doseWH,  doseWH,  doseWH,  doseFM, doseUN]

function FukushimaCost()
    # Total cost variable which will get the cost of decontaminating every pairing added.
    totalcost = 0
    for i in 2:size(input,1) # loop over rows
        for j in 1:size(input,2) # and over columns
            totalcost = totalcost + input[i,j]*solveOpt(areamapping[i-1],dosemapping[i-1],input[1,j]) 
       end
    end
    return totalcost
end
print("The total cost of decontaminating Fukushima is: \$", FukushimaCost())
The total cost of decontaminating Fukushima is: $3.2183880646568525e8

3.B. Pareto Tradeoff

Here, we investigate a tradeoff of the cost and decontamination factor of a certain area.

First, we must slightly modify our model to remove the minimum decontamination requirement and change the objective to trade off between the decontamination and cost.

In [5]:
using JuMP
function solveTradeoff(area, doseContrib, λ)
    m = Model()

    # Set 1 binary variable per each method in each surface
    @variable(m, use[1:length(n_method)], Bin)

    # Use at most one method per surface. (We can use zero if other surfaces make up for it.)
    @constraint(m, sum{use[i], i in (n_roads)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_rsoil)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_osoil)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_trees)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_walls)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_roofs)-1} <= 1)
    @constraint(m, sum{use[i], i in (n_inside)-1} <= 1)

    # Limit methods to ones the user allows
    @constraint(m, use .<= methodUse) 
                            
    # Limit methods based on time.
    for i in 1:length(n_method)
        @constraint(m, use[i]*givenTime <= methodTime[i])
    end
    
    # Want to trade off cost and decontamination
    @objective(m, Min, sum{use[i]*methodCost[i]*area[findfirst(surfaces,methodSurf[i])], i in 1:length(n_method)} -
    λ*sum{use[k]*methodDF[k]*doseContrib[findfirst(surfaces,methodSurf[k])], k in 1:length(n_method)})
    status = solve(m)
    # Sum up the costs and decontamination factors
    cost = 0
    decon = 0
    for i in 1:length(n_method)
        cost = cost + getvalue(use[i])*methodCost[i]*area[findfirst(surfaces,methodSurf[i])] 
        decon = decon + getvalue(use[i])*methodDF[i]*doseContrib[findfirst(surfaces,methodSurf[i])]
    end
    return (cost, decon, getobjectivevalue(m))
end;

Next we build a Pareto curve of our results to help visualize the tradeoff between decontamination factor and cost. Cost is much, much larger, so we use a very large tradeoff parameter $\lambda$.

In [6]:
# We'll construct a pareto curve 
points = 100
cost = zeros(points)
decon = zeros(points)
# Use an enormous log spacing, since the cost can go up really dramatically with the decontamination factor
for (j,λ) in enumerate(logspace(5,7,points))
    # We'll look at apartment complexes, as there weren't significant numbers in Fukushima
    (cost[j],decon[j],uu) = solveTradeoff(AreaAP7, doseAP, λ) 
end

# Plot the tradeoff
using PyPlot
PyPlot.svg(true)
figure(figsize=(12,4))
plot( cost, decon, "b.-")
xlabel("cost")
ylabel("decontamination factor");

4. Results and discussion

In our analysis of the radiation release at Fukushima, we see that our optimiziation gets a cost of roughly $320 million. However, let us assume that there is a shortage of small machinery able to remove topsoil, and so it is unavailable for decontamination.

In [7]:
# Can't use machinery to dig topsoil
methodUse[5] = 0
print("New cost to decontaminate Fukushima: \$", FukushimaCost())
New cost to decontaminate Fukushima: $5.842115213296392e8

By this, we can see that factors that influence the available decontamination procedures can have a large effect on the cost of decontamination. Thus, we can see that our limited selection of decontamination procedures heavily influence our eventual decontamination costs, and the assumptions made about their effectivenesses and costs.

The topic of tradeoffs is explored further in the Pareto curve visible in section 3B. Through the pareto curve we can see the dramatic increases in cost to decontaminate even just a single square kilometer of apartment buildings as the required decontamination factor increases. We can also see an interesting phenomenon where there are only a few points that the curve takes, due to the limited number of different combinations of decontamination methods.

It is worth noting that all of the costs we've seen have been on the low side compared to what the Japanese government has been spending on cleanup. This is not necessarily due to inefficiencies in the Japanese government, but instead due to modeling constraints, as this model does not consider many externalities like disposal and treatment of waste, mangerial and logistical costs, or compensation for relocated peoples. Also, there will be many other inherent constraints for a given situation that the user would need to consider that are not modelled here, which will naturally also increase costs.

5. Conclusion

We have developed an optimization routine that allows decision makers to find optimal strategies for decontamining areas for human rehabitation, and demonstrated it on the Fukushima radiation release. A major potential future direction would be to add the ability to analyze the financial considerations at a larger level, allowing a decision maker with a fixed amount of money to prioritize areas for decontamination and others for condemnation (where former residents would never return, and no decontamination measures would be necessary). Furthermore, potential health effects could be modeled based on the various decontamination methods, allowing for a more holistic approach to finding optimal decontamination strategies that does not assume all people recieve the same doses doing the same activities.