% SIMIODETemplate-TeacherORStudent.tex and simiode.cls for LaTeX2e
% Version 29 July 2015.
%
% This class (simiode.cls) is a modification of that used by
% the journal PRIMUS - Problems, Resources, and Issues in
% Mathematics Undergraduate Studies.
%
% See http://www.tandfonline.com/loi/upri .
%
% The class was developed by  Marc Frantz with suggestions of the
% PRIMUS Editorial Board, especially John Thoo.
%
% This template is to be used with the  class file "simiode.cls".
%
% The latest version of this class and template files
% can be found at
%
%       <http://www.SIMIODE.org> in the Author section.
%
% A "Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
%
\documentclass{simiode}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Basically the Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review (but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.
% Rename the pdf after you compile each final version.
% This option was created  by Eric Sullivan, Carroll College, Helena MT.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{hyperref}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding double spacing %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Manuscripts should be submitted doubled spaced.  This is
% accomplished via the \baselinestretch command below.
% However, you may obtain an idea of how your article will
% look in final form, so that you can better plan your
% presentation.  To see your manuscript formatted very
% similarly to the on-line published format of  SIMIODE,
% comment the following line out by inserting a percent sign % (%) at its beginning.

\renewcommand{\baselinestretch}{1.3}

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE regarding macros %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Any macros defined for your paper should be contained in
% the top matter. Likewise, any environment definitions such
% as \newtheorem or \newenvironment should also go in the
% top matter.
%
% Do not \input your macros as separate files within your
% final version because more files make it harder for the
% editors and users to keep track of and/or modify and
% append your paper.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Top matter is above and
% the beginning of the document below.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\newpage
% You only need to change the NAME OF SCENARIO in the next line and delete the footnote. See above for instructions on toggling from Teacher Version to Student Version.
\title{\ifnum\TeacherEdition=1 TEACHER VERSION  \else STUDENT VERSION \fi\\CONDENSATION OPTIMIZATION}




\markboth{Condensation}{Condensation}






\author{
\begin{minipage}[t]{4.5 true cm}
Chris McCarthy\\ 
Borough of Manhattan Community College\\
New York NY 10007 USA\\
\teacher{mccarthybmcc@yahoo.com}
\end{minipage}
\ \qquad \
\begin{minipage}[t]{4.5 true cm}
Brian Winkel, 
Director SIMIODE\\
Chardon OH 44024 USA\\
BrianWinkel@simiode.org
\end{minipage}
}

\makeStitlePDFLaTex


\begin{abstract}
We seek to  optimize a condensation process which is modeled by  a simulation using  the random motion of 200 particles in a 50 by 50 square in which a particle bounces off the two vertical and top walls and condenses  on the bottom wall. The model is  a differential equation in which parameters need to be estimated. We seek to determine the optimal step size so as to attain 90\% condensation level  while minimizing  costs associated with duration of the process and energy costs associated with varying the step size.
\end{abstract}


\section*{SCENARIO DESCRIPTION}

\subsection*{Prologue Describing Simulation}

Consider  the random motion of 200 particles in a 50 by 50 square with vertices $(0,0), (0,50)$, $(50,50), (0,50)$ in a plane. At each iteration, $n = 0, 1, 2, \ldots$ each particle moves a fixed step size  in one direction only, due north, due east, due south, or due west, and it must move.

Particles are initially randomly distributed on points with integer coordinates, but not on the boundaries, i.e. from within the set $S = \{(m,n)\,\vert\, m,n \in \{1, 2, 3, 
\ldots, 49\}\,\}$ and at each iteration each particle moves. The number of condensed particles is tallied at each iteration.

If a particle comes in contact with the west, east, or north wall it bounces off by returning to the position just before the bump.  If the particle comes in contact with the south or bottom wall it condenses and stays exactly at the position of contact, thus depleting the number of particles which are still randomly moving in the square. In Figures~\ref{fig:1-47-SimPlotBoth1}-\ref{fig:1-47-SimPlotBoth9000} we show a snapshot of several iterations of the particles and the plot of the number of condensed particles at each iteration using a step size of 1 unit.
We reproduce here the activities from Modeling Scenario 1-47-S-Condensation-StudentVersion for completeness.
We provide several animations for different runs of the simulation  and we offer several data sets from such simulations.

\begin{enumerate}
\item Offer up a difference or differential equation which models $y(t)$ the number of condensed particles at iteration $t$.
\item Using the data set with step size 1 in Table~\ref{table:dataset} (and shown in Figure~\ref{fig:1-47-dataset}) or from one of the data sets provided estimate the parameter(s) in your model.
 
\begin{table}[!h]
\small
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|} \hline
 Iteration & 0  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  1000  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  2000  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  3000  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  4000  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  5000  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  6000  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  7000  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  8000  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm} 
   9000  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  10000  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  11000  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  12000  \hspace*{-.2 true cm} \\ \hline
 \# Condensed & 1  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  75  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  103  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  121  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  129  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  136  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  150  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  160  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  170  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  173  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  179  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm} 
   182  \hspace*{-.2 true cm} &  \hspace*{-.2 true cm}  188  \hspace*{-.2 true cm} \\ \hline
\end{tabular}
\normalsize
\parbox{8.8 true cm}{\caption{\label{table:dataset} Sampled data  at increments of 1000 iterations from a simulation of duration 10,000 iterations, with step size 1.}}
\end{center}
\end{table}

\begin{figure}[!h]
\begin{center}
\includegraphics[width=9cm]{1-47-dataset.eps}


\parbox{9 true cm}{\caption{\label{fig:1-47-dataset} Plot of data in Table~\ref{table:dataset} of number of condensed particles at iterations, 1000, 2000, 3000, \ldots, 11,000, 12000.}}
\end{center}
\end{figure}

\item Confirm in some way the correctness of your model.
\end{enumerate}


\begin{figure}
\centering
\hspace{-1.5 true cm}\includegraphics[width=16cm]{1-47-SimPlotBoth1.eps}
\vspace*{-.15 true cm}
\parbox{13 true cm}{\caption{\label{fig:1-47-SimPlotBoth1} On the left we see the particles in the box at iteration 1 with 1 of the particles condensed  on the bottom wall and on the right we see a plot of the accumulated number of condensed particles at given iteration.}}


\vspace*{.3 true cm}
\centering
\hspace{-1.5 true cm}\includegraphics[width=16cm]{1-47-SimPlotBoth6000.eps}
\vspace*{-.15 true cm}
\parbox{13 true cm}{\caption{\label{fig:1-47-SimPlotBoth6000} On the left we see the particles in the box at iteration 6000 with 161 of the particles condensed  on the bottom wall and on the right we see a plot of the accumulated number of condensed particles at given iteration.}}


\vspace*{.3 true cm}
\centering
\hspace{-1.5 true cm}\includegraphics[width=16cm]{1-47-SimPlotBoth9000.eps}
\vspace*{-.15 true cm}
\parbox{13 true cm}{\caption{\label{fig:1-47-SimPlotBoth9000} On the left we see the particles in the box at iteration 9000 with 187 of the particles condensed  on the bottom wall and on the right we see a plot of the accumulated number of condensed particles at given iteration.}}
\end{figure}

\newpage
\subsection*{Optimization}

In industrial applications it is often desirable to extract or separate particles, e.g., purification or collection for sale or re-use. Indeed, our condensation efforts are but one example of extraction.

Cost is a significant concern in all industrial efforts and we would like to examine our process in light of minimizing costs associated with condensation of particles.  Accordingly, we can alter variables for a goal of separating out 90\% of our 200 particles, i.e. 180 particles.

A reasonable assumption is that our process of condensation needs more heat for the particles to  take a larger step size (more energy means more excitation) and so doing the condensation will be faster, i.e. the time for condensation will be reduced.

For our process let us ignore the ``purchase costs" of particles and possible 	``reimbursements" for unused particles, the 20 particles not condensed in each operation. Also we will ignore set up and tear down costs for the apparatus.

Suppose the energy cost to produce a step size, $h$, in our particle motion is ${\rm EnergyCost}(h)$. A reasonable relationship between energy and step size might be linear, i.e. ${\rm EnergyCost}(h) = \alpha h$, meaning for each additional unit of increase in step size our costs go up by $\alpha$.

Further, suppose the cost to run the experiment for time ${\rm T}(h)$ as a function of $h$, step size, i.e. the time until, say, 90\% or 180 of our particles are condensed, is ${\rm TimeCost}(h)$.

Then the total cost to run our process could be written as a function of our step size, $h$,

\begin{equation}\label{eq:totl}
				{\rm TotalCost}(h) = {\rm EnergyCost}(h) + {\rm TimeCost}(h)\, .
\end{equation}

We will need to determine ${\rm T}(h)$, the time our process takes to condense 180 of our 200 particles as a function of step size, $h$, and then   ${\rm TimeCost}(h) = \beta {\rm T}(h)$ and ${\rm EnergyCost}(h) = \alpha h$  where $\beta$ is a constant for the  cost per unit of time of running our process and $\alpha$ is the incremental cost in energy in order  to raise the step size one unit.

Suppose our cost constants are  $\beta = 0.0012$ and $\alpha = 1.5$.



\subsubsection*{Modeling ${\rm T}(h)$} 

We ran the condensation simulation using several different step sizes, $h = 1, 2, 3, 4, 5, 6$. For each step size we ran 50 simulations and we averaged the time it took for a simulation to reach exactly 180 condensed particles. So, for example, with step size $h = 3$ the average amount of time it took the simulation to reach exactly 180 condensed particles was 955.34 units of time. Complete data is offered in Table~\ref{table:Tdata}.

\begin{table}[!h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline
 Step size ($h$) & 1 & 2 & 3 & 4 & 5 & 6 \\ \hline
 Duration of simulation (${\rm T}(h)$) & 8245.34 & 2161.48 & 955.34 & 551.68 & 371.57 & 264.24 \\ \hline
\end{tabular}
\parbox{12 true cm}{\caption{\label{table:Tdata}  Empirical data on step size $h$ and ${\rm T}(h)$, the average time (of 50 runs for each step size) for the condensation process to yield exactly 90\%  of the original 200 particles or 180 particles condensed.}}
\end{center}
\end{table}

\subsubsection*{Incidental - Modeling the relationship between step size $h$ and parameter $k$} 

We used the   differential equation model, (\ref{eq:DEModel}), quite successfully, to model, $y(t)$, the number of condensed particles at time $t$

\begin{equation}\label{eq:DEModel}
y'(t) = k (200 - y(t)\, , \quad y(0) = 0\, .
\end{equation}


Data on step size and the resulting parameter $k$ using model (\ref{eq:DEModel}) was collected from 10 runs of the simulation for condensation, each of duration 10,000 iterations, to determine 10 values of $k$. These values were averaged for each step size and the results are summarized in Table~\ref{table:kVal}. 

\begin{table}[!h]
\small
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}\hline 
 Step size ($h$)  \hspace{-.2 true cm} & \hspace{-.2 true cm}   1  \hspace{-.2 true cm} & \hspace{-.2 true cm}   2  \hspace{-.2 true cm} & \hspace{-.2 true cm}   1.5  \hspace{-.2 true cm} & \hspace{-.2 true cm}   2.5  \hspace{-.2 true cm} & \hspace{-.2 true cm}   3  \hspace{-.2 true cm}  \\ \hline
 Parameter $k$ \hspace{-.2 true cm} & \hspace{-.2 true cm}   0.0003491  \hspace{-.2 true cm} & \hspace{-.2 true cm}   0.001241  \hspace{-.2 true cm} & \hspace{-.2 true cm}   0.0006743  \hspace{-.2 true cm} & \hspace{-.2 true cm}   0.001699  \hspace{-.2 true cm} & \hspace{-.2 true cm} 0.00265 \hspace{-.2 true cm}\\ \hline
\end{tabular}

\vspace*{.3 true cm}
\hspace*{1 true cm}\begin{tabular}{|c|c|c|c|c|c|c|}\hline 
\hspace{-.2 true cm}   3.5  \hspace{-.2 true cm} & \hspace{-.2 true cm}   4  \hspace{-.2 true cm} & \hspace{-.2 true cm}   4.5  \hspace{-.2 true cm} & \hspace{-.2 true cm}   5  \hspace{-.2 true cm} & \hspace{-.2 true cm}   5.5  \hspace{-.2 true cm} & \hspace{-.2 true cm}   6  \hspace{-.2 true cm}\\ \hline
 \hspace{-.2 true cm}   0.003465
    \hspace{-.2 true cm} & \hspace{-.2 true cm}   0.004958  \hspace{-.2 true cm} & \hspace{-.2 true cm}   0.00605  \hspace{-.2 true cm} & \hspace{-.2 true cm}   0.006953  \hspace{-.2 true cm} & \hspace{-.2 true cm}   0.008286  \hspace{-.2 true cm} & \hspace{-.2 true cm}   0.009922  \hspace{-.2 true cm}\\ \hline
\end{tabular}



\parbox{10.5 true cm}{\caption{\label{table:kVal} Empirical  data on step size $h$ and the parameter $k$ from model (\ref{eq:DEModel}) which resulted from 10 simulations of duration 10,000 iterations for each step size. For each step size the 10 obtained values of $k$ were averaged.}}
\end{center}
\end{table}


\subsection*{Questions to Address}

\begin{enumerate}
\item From the data in Table~\ref{table:Tdata} determine a functional relationship for ${\rm T}(h)$, the average time (of 50 runs for each step size) for the condensation process to yield exactly 90\%  of the original 200 particles or 180 particles condensed.
\item From the data in Table~\ref{table:kVal} determine a functional relationship  between step size $h$ and the corresponding parameter $k$. Explain the reasonableness of this relationship.
\item  Construct the  Total Cost function ${\rm TotalCost}(h)$ from (\ref{eq:totl}).
\item Determine the step size, $h$, which gives minimum ${\rm TotalCost}(h)$.

\end{enumerate}

\teacher{\newpage
\section*{NOTES FOR TEACHERS}

Let us address  question (1) ``From the data in Table~\ref{table:Tdata} determine a functional relationship for ${\rm T}(h)$, the average time (of 50 runs for each step size) for the condensation process to yield exactly 90\%  of the original 200 particles or 180 particles condensed.''

When we attempt to fit a function to the data we find a very nice empirical fit with the function 

\begin{equation}\label{eq:funcfit}
{\rm T}(h) = \frac{8277.08}{h^2}\, ,
\end{equation}
and this can be seen in the plot offered in Figure~\ref{fig:1-47a-StepTime}.

\begin{figure}[!h]
\begin{center}
\includegraphics[width=12cm]{1-47a-StepTime.eps}


\parbox{11.5 true cm}{\caption{\label{fig:1-47a-StepTime} Plot of data for and best fit  relationship for ${\rm T}(h)$, the average time (of 50 runs for each step size) for the condensation process to yield exactly 90\%  of the original 200 particles or 180 particles condensed.
}}
\end{center}
\end{figure}


\subsection*{Addressing  the relationship between step size $h$ and parameter $k$}

Let us consider the incidental question of modeling the relationship between step size $h$ and parameter $k$.  We see the result of attempting three different models: (1) (red)  $k = a h + b h^2$, (2) (blue) $k = A e^{B h}$, and (3) (green) $k = A h e^{B h}$ in Figure~\ref{fig:1-47a-StepParameter}.

\begin{figure}[!h]
\begin{center}
\includegraphics[width=12cm]{1-47a-StepParameter.eps}


\parbox{11.5 true cm}{\caption{\label{fig:1-47a-StepParameter} Plot of data from Table~\ref{table:kVal} used to determine a functional relationship  between step size $h$ and the corresponding parameter $k$. Three models were compared: (1) (red)  $k = a h + b h^2$, (2) (blue) $k = A e^{B h}$, and (3) (green) $k = A x e^{B h}$.
}}
\end{center}
\end{figure}



From a summary of our  analysis in Table~\ref{table:compare}, it would appear that the quadratic model, $k = a h + b h^2$,  has smallest least sum of square error value and so would appear to be best model for parameter value or rate constant $k$ vs step size.


\begin{table}[!h]
\begin{center}
\begin{tabular}{|c|c|} \hline
Model & SSE \\ \hline
$k = a h + b h^2$ & $3.19484\times 10^{-7}$ \\ \hline
$k = A e^{B h}$ & $2.91967\times 10^{-6}$ \\ \hline
$k = A x e^{B h}$ & $8.79646\times 10^{-7}$ \\ \hline
\end{tabular}

\parbox{8 true cm}{\caption{\label{table:compare} A comparison of the least sum of square errors (SSE) for three models used to determine a functional relationship  between step size $h$ and parameter $k$.}}
\end{center}

\vspace*{-.5 true cm}
\end{table}


It should be noted that if we take step sizes, say of $h = 0$, then in $y'(t) = k (200 - y(t))$ we would expect $k = 0$, for there would be no movement and hence no increase at all in the number of condensed particles. That means the quadratic $k = a h + b h^2$ would be better than the pure exponential $k = A e^{B h}$ as it passes through the origin $(0,0$), but we could use  $k = A h e^{B h}$ which also gives $k = 0$ for $h = 0$. However, the SSE is smallest for the quadratic form, $k = a h + b h^2$. Thus we see that as step size, $h$, increases the rate constant $k$ increases by the square of the step size essentially.  This means a doubling of the step size quadruples the speed of the condensation, something to consider in keeping costs down.



\subsubsection*{Cost Optimization}

We now turn to the total cost of the operation as a function of step size $h$ as we have all the pieces for  ${\rm TotalCost}(h) = {\rm EnergyCost}(h) + {\rm TimeCost}(h)$, namely: ${\rm TimeCost}(h) = \beta {\rm T}(h)$ and ${\rm EnergyCost}(h) = \alpha h$ where ${\rm T}(h) =   \frac{8277.08}{h^2}$. Let use $\gamma$  for our value of   8277.08  to be more general.

\begin{equation}\label{eq:speccosts}
{\rm TotalCost}(h) =  {\rm EnergyCost}(h) + {\rm TimeCost}(h) =  \alpha h+ \beta {\rm T}(h) = \alpha h+ \beta \frac{\gamma}{h^2}
\end{equation}

Now to minimize ${\rm TotalCost}(h)$ we take the first derivative of ${\rm TotalCost}(h)$ with respect to $h$ and set it equal to zero, whereupon we solve for our optimal value of $h$. This results in 

\begin{equation}\label{eq:opth}
h_{\rm opt} = \left(\frac{2 \cdot \beta \cdot \gamma}{\alpha}\right)^{1/3}\, .
\end{equation}

And when we substitute our optimal  value of $h_{\rm opt}$ into the total cost ${\rm TotalCost}(h)$ equation (\ref{eq:speccosts}) we find:

\begin{equation}\label{eq:opth}
{\rm TimeCost}(h_{\rm opt}) = 3 \left(\frac{\alpha^2 \cdot \beta \cdot \gamma}{4}\right)^{1/3}\, .
\end{equation}

So, finally putting in our numbers for $\alpha = 1.5$, $\beta = 0.0012$ and $\gamma = 8277.08$, we have a minimum cost of $${\rm EnergyCost}(h_{\rm opt}) = 5.32\, ,$$ with separate running time  cost of 
$${\rm EnergyCost}(h_{\rm opt}) = \beta {\rm T}(h_{\rm opt}) =  0.0012 \left(\frac{8277.08}{h_{\rm opt}^2}\right) = 1.77$$ and energy cost of $${\rm EnergyCost}(h_{\rm opt}) = \alpha h_{\rm opt} = 1.5 h_{\rm opt} = 3.55 \, .$$
The run time will be ${\rm T}(h_{\rm opt}) = 1479$ units of time.



We see that $${\rm TimeCost}(h_{\rm opt}) = 3 \left(\frac{\alpha^2 \cdot \beta \cdot \gamma}{4}\right)^{1/3}$$ is reasonable, for if we increase $\alpha$, our cost of energy, or $\beta$, our cost due to time, the overall cost goes up. Interestingly enough, if we multiplied $\alpha$ by a factor of 8 our costs will  quadruple, whereas multiplying $\beta$ by a factor of 8 will only double our costs.  

In summary, this analysis for our situation suggests we should set our step size at 2.37. The run time will be 1479 units of time, while the cost of running the process for that long will be 1.77. The energy cost will be 3.55 and the total cost will be 5.32.

}




\end{document}


