Report for Homework 7 - Stochastic Modeling

Write Python code to model the following:

Water Treatment System Design

You want to design a simple plant-based water treatment plant for removing a nutrient from a wastewater stream. The system contains a fixed amount of plant material, taking up nutrients according to first-order reaction kinetics (i.e. the rate of uptake is a linear function of N. Assuming the system is system is homogeneous in space (i.e. a completely mixed reactor, the governing equation for this system is:

dNdt=fV(NinN)NdNdt=fV(NinN)N

where:

N = nutrient concentration (g/m3 = mg/L)

f = flowrate (m3/day)

V = reactor volume (m3)

Nin = N concentration in Influent (g/m3)

= first-order rate constant (m3/g-day)

For this analysis, we will assume that the 1) flow rate, and 2) influent nutrient concentration Nin are variable. Specifically, we will assume that flow rate is distributed according to a lognormal distribution, and Nin is distributed according to a normal distribution. We also have observations of flow rates and Ninvalues we will use to parameterize these distributions.

We will only concern ourselves with steady state results from our model for this problem.

Model Parameters

Model parameter values are:

Parameter	Value	Units
alpha	0.004	m3/g-day
V	500	m3

Your Tasks:

Derive the model for steady-state. Write a function that takes as input arguments a flowrate and an influent N concentration and returns the resulting steady state value of N. Report the steady state N concentration for a flow rate of 4.0 m3/day, NIn=10.0 mg/L.

Read in the flow rate and influent concentation data and fit a probability distribution to the data. Use lognorm for flow rate, and norm for the influent N concentration. Report the distribution parameters for each distribution. Generate a histogram, overlayed with the fitted PDF, for each variable.

Use your steady-state model to conduct a Monte-Carlo simulation:

Generate 100 random flow values and 100 random Nin value by sampling from the distributions developed above.

Solve the steady state model for the 100x100 flow rate and Nin combinations, collecting the steady state N values and corresponding removal efficiency for each combination of flow rate and Nin, placing each result in an array

Fit the resulting steady state N's to a weibull_max distribution.

Show the generated results in a histogram on a plot, along with the PDF of the distribution.

Turn in: Turn in a Jupyter notebook with the the following:

A function implementing the steady-state solution to the model.

The distribution parameters for each of the observed variables in Task 2 above. In addition, generate two plots showing the results of you analysis for each variable. Each plot should show 1) a histogram of the observated data, and 2) the pdf of thefitted distribution function. The first plot should report results for flowrate, the second for Nin.

The distributions of steady-state N and Removal Efficiencies resulting from your Monte-Carlo analysis, depicted similarly to your plots in (2) above.

Result from answering the five questions posed below in the last markdown block:

What is the likelihood of observing flows > 15 m3/day?

What is the likelihood of observing influent concentration between 9 and 11 mg/L?

How often would removal efficiencies be expected to drop below 50 percent?

Your design goal is to maintain effluent N concentrations below 8.0 mg/L. How often would you expect this condition to be satisfied?

What is the minimum effluent N concentration could I expect to achieve 95 percent of the time?

Implement the Steady State Model

First, derive the steady state solution to the model. Then, implement the steady state model in a function GetSteadyStateN(_flowRate, _NInf) that returns the steady state value of N and removal efficiency for the specified flowrate and influent N concentration.

Test your GetSteadyStateN() with the following inputs: flowRate=4.0 m3/day, NIn=10.0 mg/L. It should produce a steady state N value of 6.67 mg/L and a removal efficiency of 33 percent. Output your result for this combination of inputs in your code below.

In [ ]:

import numpy as np

Develop Probability Distributions for Flowrate, Influent N Concentrations

Access the dataset at http://explorer.bee.oregonstate.edu/Topic/Modeling/Data/WTParams.csv - it contains two columns, 'FlowRate' and 'Influent N', each of which contains multiple observations of these two variables. This will be the data we use to parameterize the two distributions.

Hint: to generate a histogram, use the matplotlib.pyplot function hist, setting the normed keyword argument to True, e.g.

plt.hist(Y, bins=25, normed=True, color='g', label='Histogram') # 'Y' is the data to plot

Report the distribution parameters for each distribution. Additionally, generate a figure with two subplots: the left plot should show the histogram and PDF for the flow rate data, the right the histogram and PDF for the Nin data.

In [ ]:

import numpy as np

import pandas as pd

from scipy.stats import lognorm

from scipy.stats import norm

import matplotlib.pyplot as plt

Run a Monte-Carlo Simulation

Generate 100 random flow values and 100 random Nin values by sampling from the distributions developed above.

Solve the steady state model y calling your GetSteadyStateN() function, using the 100x100 flow rate and Nin combinations, calculating steady state removal efficiencies from your steady-state model results. To do this, use a nested loop - the outer looping through your random flowrate value, and the inner through your random Nin values. Store the steady state Removal Efficiencies for each model run in an array (it will end up with 100x100=10000 of these values at the end of your simulation.

Fit the resulting steady state Removal Efficiencies to a weibull_max distribution.

Generate a figure with two subplots, side-by-side. The left panel should include a histogram overlaid with your PDF. The right panel should show the CDF. Use a range of [0,1] for your x-axis. Label appropriately

In [ ]:

import numpy as np

from scipy.stats import lognorm

from scipy.stats import norm

from scipy.stats import weibull_max

import matplotlib.pyplot as plt

Answer the following questions:

Based on the results obtained in your Monte-Carlo simulation, determine:

What is the likelihood of observing flows > 15 m3/day?

What is the likelihood of observing influent concentration between 9 and 11 mg/L?

How often would removal efficiencies be expected to drop below 50 percent?

Your design goal is to maintain effluent N concentrations below 8.0 mg/L. How often would you expect this condition to be satisfied?

What is the minimum effluent N concentration could I expect to achieve 95 percent of the time?