Cortix 2019 03Aug2019

City Criminal Justice Dynamics Example

  • This is part of the Cortix-on-Jupyter-Notebook examples.
  • You must be in a Jupyter Notebook server to run this notebook.
  • Select each of the cells below and run them sequentially (use the run button, >| on the tool bar or use the Cell option on the menu bar).
  • Alternatively, on the menu bar run all cells: Cell -> Run All.

$ \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Smtrx}{\boldsymbol{\mathsf{S}}} \newcommand{\xvec}{\boldsymbol{\mathsf{x}}} \newcommand{\vvar}{\boldsymbol{v}} \newcommand{\fvar}{\boldsymbol{f}} \newcommand{\Power}{\mathcal{P}} \newcommand{\bm}[1]{{\boldsymbol{#1}}} $

Introduction

This notebook uses Cortix to simulate the startup and shutdown of a boiling water reactor using a point-reactor model, This work is in progress, authors: Prof. Valmor F. de Almeida, UMass Lowell, Dept. of Chemical Engineering; And Austin Rotker. Cortix Group.

BWR Model

The following equations are developed for use in the point-reactor BWR model.

One-Group Energy Neutron Balance

The space-invariant neutron balance for the point-reactor model is

\begin{equation*} \frac{\text{d}n}{\text{d}t} = \frac{\rho(t)-\beta}{\ell}\, n + \sum\limits_{i=1}^{6} \lambda_i\,c_i + q(t), \end{equation*}

where the first term on the right side of the equation represents the net production of neutrons not accounting for delayed neutrons, and the second term accounts for the source of delayed neutrons considering 6 groups of delayed neutron emitters resulting from the fission of $^{235}$U nuclei. Therefore a balance of neutron emitter species is also necessary

\begin{equation*} \frac{\text{d}c_i}{\text{d}t} = \frac{\beta_i}{\ell}\, n - \lambda_i\,c_i , \ \ \ \ \ \forall \ \ \ \ \ i=1,\ldots,6. \end{equation*}

where the first term on the right side of the equation is the source of emitters as a function of the neutron number density $n(t)$, and the second term is the consumption rate of the emitter by radioactive decay obtained as a function of the product of the number concentration of the emmiter, $c_i(t)$, multiplied by its decay constant $\lambda_i$. Here the number concentration of of the $i$th emitter, $c_i$ is considered in terms of delayed neutron number density, hence the units are the same as $n$.

The current model considers the temperature reactivity coefficient, $\alpha_{T_n}$, that is,

\begin{equation*} \rho(t) = \rho_0 + \alpha_{T_n}(T_f)\,\bigl(T_f(t) - T^{(0)}\bigr), \end{equation*}

A heat balance in the static fuel dispersed into the moderator gives

\begin{equation*} \frac{\text{d}T_f}{\text{d}t} = - \frac{1}{\rho_f\,c_{pf}}\biggl(q^{'''}(t) - \frac{\dot{Q}_\text{f}}{V_f} \biggr) \end{equation*}

where the heat rate sink in the fuel/moderator is

\begin{equation*} {\dot{Q}_\text{f}} = -UA\bigl(T_f-T_c\bigr), \end{equation*}

and the nuclear heating power is given by

\begin{equation*} q^{'''}(t) = G_f\frac{\sqrt{\pi}}{2}\,\sigma_\text{fo}\,\sqrt{\frac{T^\text{(o)}}{T}}\,\biggl(\frac{\epsilon\,w_\text{fs}\,\rho_\text{fm}}{M_\text{fn}}\,i\,N_A\biggr)\,n(t)\,v_\text{(o)} . \end{equation*}

A heat balance in the flowing coolant fluid produces

\begin{equation*} \frac{\text{d}T_c}{\text{d}t} = -\frac{1}{\tau_c}\biggl(T_c-T_{in}\biggr) - \frac{1}{\rho_c\,c_{pc}}\biggl( \frac{\dot{Q}_\text{f}}{V_c} \biggr) \end{equation*}

where the heat rate source is the negative of the heat rate sink in the fuel/moderator.

Neutron density and delayed neutron precursor concentrations at steady-state are related as follows:

\begin{align*} n_\text{ss} &= -\frac{q_\text{ss}\,\ell }{\rho_\text{ss} } \\ c_{i_\text{ss}} &= \frac{\beta_i}{\lambda_i\,\ell}\, n_\text{ss} \ \ \ \ \forall \ \ \ \ i=1,\ldots,6 . \end{align*}

Moderator Temperature Feedback

The moderator temperature reactivity coefficient is the major contributor to the power coefficient in a boiling water reactor. It can be calculated by the following equation:

\begin{equation*} \alpha_\text{m} = \frac{dK}{dT_\text{m}}\frac{1}{K}\ \end{equation*}

Where $\alpha_\text{m}$ is the moderator temperature feedback coefficient, K is the neutron multiplecation factor for a finite medium, and $T_\text{m}$ is the temperature of the moderator. K is calculated as follows:

\begin{equation*} K = \frac{\eta \epsilon P F}{1 + (\tau + L^2)B^2 + \tau L^2 B^4} \end{equation*}

Plugging this in to the original moderator temperature reactivity coefficient equation leaves the following differential equation:

\begin{equation*} \alpha_\text{m} = \frac{d}{dT_\text{m}}\frac{1 + (\tau + L^2)B^2 + \tau L^2 B^4}{\eta \epsilon P F}\ \end{equation*}

With $\tau$ being the fermi age or slowing down area, $L$ being the neutron diffusion length, $B$ being buckling, $\eta$ being neutron yield, $\epsilon$ being fast fission factor, $P$ being resonance escape probability and $F$ being thermal utilization. Only $\tau$, $P$ and $F$ vary significantly with moderator temperature, with $L$ and $B$ being assumed constants over the range of temperatures.

$\tau$ is found as follows:

\begin{equation*} \tau = \frac{D_\text{f}}{\Sigma_\text{1}} \end{equation*}

with $D_\text{f}$ being the diffusion coefficient and equal to:

\begin{equation*} D_\text{f} = \frac{1}{3 \Sigma_\text{s}(1 - \mu_\text{0})} \end{equation*}

With $\Sigma_\text{s}$ being the macroscopic scattering cross section, $\sigma_\text{s} N_\text{m}$, and $\mu_\text{0}$ being a constant dependent on the material used as the moderator.

$\Sigma_\text{1}$ is the neutron removal cross section, or the probability per unit length that a neutron becomes thermalized. It is calculated as:

\begin{equation*} \Sigma_\text{1} = \frac{\xi \Sigma_\text{s}}{ln \frac{E_\text{0}}{E_\text{th}}} \end{equation*}

With $\xi$ being the average logarithmic energy decrement for a collision involving a neutron and the moderator (in the case of a light water reactor, $\xi$ = 1), $E_\text{0}$ being the energy of a neutron produced during fission (assumed to be constant at around 2 MeV), and $E_\text{th}$ being the energy of the medium that the neutron is diffusing through, which is dependent on the temperature of the moderator by the following equation:

\begin{equation*} E_\text{th} = (T_\text{m})\frac{0.0862 \ MeV}{K} \end{equation*}

With $T_\text{m}$ being in kelvin.

Assuming that only $N_text{m}$ and $E_\text{th}$ vary with the temperature of the moderator, then the differential of $\tau$ with respect to $T_\text{m}$ may be calculated as follows:

\begin{equation*} \frac{d \tau}{d T_\text{m}} = \frac{d}{d T_\text{m}} \frac{\frac{1}{3 \sigma_\text{s} N_\text{m}(1 - \mu_\text{0})}}{\frac{\xi \sigma_\text{s} N_\text{m}}{ln \frac{E_\text{0}}{E_\text{th}}}} \end{equation*}\begin{equation*} \frac{d \tau}{d T_\text{m}} = \frac{((0.0862\frac{E_\text{th}}{E_\text{0}})(3 N_\text{m} )-(ln \frac{E_\text{0}}{E_\text{th}})(6 \frac{d N_\text{m}}{d T_\text{m}})}{(3 N_\text{m})^2 (1 - \mu_\text{0})} \end{equation*}

Next, the diffusion area in the moderator before absorption, $L^2$, may be evaluated. It can be calculated through the following relation:

\begin{equation*} L^2 = \frac{1}{3 \Sigma_\text{s} \Sigma_\text{a} (1 - \mu_\text{0})} \end{equation*}

And therefore:

\begin{equation*} L = \sqrt{\frac{1}{3 \Sigma_\text{s} \Sigma_\text{a} (1 - \mu_\text{0})}} \end{equation*}

Note that L is referred to as the diffsion length whereas L^2 is the diffusion area.

Taking the derivative of the diffusion length, L, with respect to time, and noting that only the number density of the moderator varies with time, and that the macroscopic cross sections $\Sigma$ can be rewritten as $\sigma N_\text{m}$, we obtain the differential of diffusion length with respect to time:

\begin{equation*} \frac{dL}{dT_\text{m}} = \frac{1}{2 \sqrt{\frac{-2 \frac{dN_\text{m}}{dT_\text{m}}}{3 \sigma_\text{s} \sigma_\text{a} N_\text{m}^3 (1 - \mu_\text{0})}}} \end{equation*}

Next, the terms in the four factor formula, $\eta \ , \ \epsilon \ , \ P, \ F$ must be developed and the derivatives of the ones which vary with moderator temperature must be calculated.

$\eta$ and $\epsilon$ do not vary with moderator temperature, and are assumed to have constant values of 2.02 and 1.03.

$P$, the resonance escape, is calculated via the following equation:

\begin{equation*} P = e^{- \frac{N_\text{f} V_\text{F} I}{\xi \Sigma_\text{s M} V_\text{M}}} \end{equation*}

In the above equation, $N_\text{f}$ is the number density of the fuel, $V_\text{F}$ is the volume of the fuel, $I$ is the resonance escape integral, and $V_\text{M}$ is the volume occupied by the moderator.

With only the macroscopic scattering cross section of the moderator varrying with temperature of the moderator (volume of coolant in the core is assumed constant, only the density changes), the derivative of the resonance escape probability with respect to the temperature of the moderator can be calculated as:

\begin{equation*} \frac{dP}{dT_\text{m}} = \frac{-N_\text{f} V_\text{F} I \sigma_\text{s} \frac{dN_\text{m}}{dT_\text{m}}}{(\xi \Sigma_\text{s M} V_\text{M})^2}e^{- \frac{N_\text{f} V_\text{F} I}{\xi \Sigma_\text{s M} V_\text{M}}} \end{equation*}

The thermal utilization, $F$, is calculated as follows:

\begin{equation*} F = \frac{\Sigma_\text{aF}}{\Sigma_\text{aF} + \Sigma_\text{aM}} \end{equation*}

With $\Sigma_\text{aF}$ being the macroscopic absorption cross section of the fuel and $\Sigma_\text{aM}$ being the macroscopic absorption cross section of the moderator. As only the macroscopic absorption cross section of the moderator will vary with the temperature of the moderator, the derivative of thermal utilization with respect to the temperature of the moderator can be written as follows:

\begin{equation*} \frac{dF}{dT_\text{m}} = \frac{-(\Sigma_\text{aF})(\sigma_\text{aM} \frac{dN_\text{m}}{dT_\text{m}})}{(\Sigma_\text{aF} + \Sigma_\text{aM})^2} \end{equation*}

The final preliminary differential equation is the one that must be developed for the differential of the moderator number density with respect to the temperature of the moderator. The change in number density of the moderator with respect to temperature will be proportional to the change in the mass density of the moderator with respect to temperature, by the molar mass of the moderator (in this case, 18 g/mol) and avogadro's number. The derivative of the mass density with respect to the temperature of the moderator of water is a function of a forty term polynomial used by the IAPWS 97 standards to approximate the behavior of water at different temperatures and pressures, that is:

\begin{equation*} \frac{d \rho}{dT_\text{m}} = \frac{d}{dT_\text{m}}(IAPWS(T)) \end{equation*}

The derivative can be approximated by using scipy's inbuilt differentiation function. From there, the differential of mass density with respect to moderator temperature may be related to the differential of number density with respect to moderator temperature by the molar mass and avogadro's number:

\begin{equation*} \frac{dN_\text{m}}{dT_\text{m}} = (18 g/mol)(N_\text{A})(\frac{d \rho}{dT_\text{m}}) \end{equation*}

Now that all previous relationships have been developed, the moderator temperature feedback coefficient can be calculated. Recall the original equation for the moderator temperature feedback coefficient:

\begin{equation*} \alpha_\text{m} = \frac{d}{dT_\text{m}}\frac{1 + (\tau + L^2)B^2 + \tau L^2 B^4}{\eta \epsilon P F}\ \end{equation*}

Taking the derivatives with respect to moderator temperature yields the following differential:

\begin{equation*} \alpha_\text{m} = \frac{(\eta \epsilon P F)((\frac{d \tau}{dT_\text{m}})(B^2 + L^2 B^4)+(\frac{dL}{dT_\text{m}})(2L B^2 + 2L \tau B^4))-(1 + (\tau + L^2)B^2 + \tau L^2 B^4)(\eta \epsilon)(F \frac{dP}{dT_\text{m}} + P \frac{dF}{dT_\text{m}})} {(\eta \epsilon P F)^2} \end{equation*}

Turbine model

Temperature and pressure of the coolant leaving the reactor is given by temp_c and pressure; these two variables are used with the IAPWS package (under the IAPWS-97 standards) to determine the specific enthalpy of the steam leaving the reactor. The turbine expands any steam passing through it to 0.75 kPa. This process is NOT assumed to be isentropic, and instead operates at 80% of maximum isentropic efficiency.

We start with a basic heat balance on the turbine:

\begin{equation*} Q + W_\text{s} = \Delta H + \Delta K + \Delta Z \end{equation*}

With Q being the heat extracted from the working fluid, W being the work done by the working fluid (and hence being positive if the turbine is producing useful work), and with $\Delta$H, $\Delta$K and $\Delta$Z representing the changes in enthalpy, kinetic energy and gravitational potential energy of the fluid, respectively. It is assumed that the turbine operates adiabatically and with no change in the kinetic or gravitational potential energy of the fluid as it moves through the turbine. Hence, the equation reduces to:

\begin{equation*} W_\text{s} = \Delta H \end{equation*}

If the turbine is producing useful work, then W is positive, and hence $\Delta$H must be positive as well. The enthalpy of the fluid will decrease as it passes through the turbine and is expanded to a lower temperature and pressure, hence the equation can be rewritten as:

\begin{equation*} W_\text{s} = \Delta H = H_\text{inlet} - H_\text{outlet} \end{equation*}

The inlet steam is at the same state of the steam leaving the top of the reactor, I.E. it has a pressure $P_\text{reactor}$ and the same temperature of the coolant, $T_\text{coolant}$. The turbine expands whatever is taken into the inlet to a final pressure of 0.75 kPa and a temperature $T_\text{runoff}$ which varies based on the properties of the steam that was fed to the turbine. Hence, the work balance equation can be re-written as:

\begin{equation*} W_\text{s} = H_\text{inlet}(T_\text{coolant}, P_\text{reactor}) - H_\text{outlet}(T_\text{runoff}, 0.75 kPa) \end{equation*}

Ideally, final temperature, $T_\text{runoff}$, would be the bubble point temperature of the runoff at 0.75 kPa (~ 3 C$^{\circ}$). However, the turbine does not operate isentropically, and instead is assumed to operate at 80% of its maximum isentropic efficiency.

We denote the maximum possible isentropic work of the turbine, $W_\text{isentropic}$, as the amount of work that would be done if the steam entering the turbine was expanded to a pressure of 0.75 kPa and a temperature of 3 C$^{\circ}$.

\begin{equation*} W_\text{isentropic} = H_\text{steam}(T_\text{coolant}, P_\text{reactor}) - H_\text{water}(0.75 kPa, 3 ^oC) \end{equation*}

The efficiency of the turbine is then given in terms of % of maximum isentropic efficiency. For example, with an isentropic efficiency of 80%:

\begin{equation*} W_\text{real} = 0.8 \, W_\text{isentropic} \end{equation*}

And, by an energy balance, we can calculate the ending enthalpy of the liquid runoff at 0.75 kPa, $H_\text{real liquid}$:

\begin{equation*} W_\text{real} = H_\text{steam}(T_\text{coolant}, P_\text{reactor}) - H_\text{real liquid}(T_\text{runoff}, 0.75 kPa) \end{equation*}

With the ending enthalpy of the liquid runoff known, and its pressure fixed by the output pressure of the turbine (0.75 kPa), the runoff temperature of the liquid from the turbine can be calculated using the IAPWS97 backwards solver functions.

Condenser Model

The condenser takes the runoff from the turbine and condenses it to a subcooled liquid at the current operating pressure of the reactor. The work required to do this is calculated as follows:

\begin{equation*} W_\text{condenser} = \frac{H_\text{turbine runoff}(T_\text{runoff}, 0.0035 kPa) - H_\text{subcooled liquid out}(T_\text{subcooled}, P_\text{reactor})}{efficiency}) \end{equation*}

Where the enthalpy of the subcooled liquid out is given by:

\begin{equation*} H_\text{subcooled liquid out}(T_\text{subcooled}, P_\text{reactor}) = H_\text{saturated liquid}(T_\text{saturated}, P_\text{reactor}) * (subcooling) \end{equation*}

With subcooling being the % by which the enthalpy of the subcooled liquid is lowered from the enthalpy it would have if it were at its boiling point. This % is constant, and is based on the % of subcooling that exists at steady state with a reactor pressure of 7 mPa and a condenser runoff temperature of 220 degrees celcius:

\begin{equation*} subcooling = (1 -\frac{H(T_{saturated}, 7 mPa, x=0) - H(220 C, 7 mPa}{H(T_{saturated}, 7 mPa, x=0}) \end{equation*}

Write the run context

In [1]:
# Import various packages; must have the Cortix repository installed

import matplotlib.pyplot as plt
import scipy.constants as const

# Leave this block here for Azure
try:
    import cortix
except ImportError:
    print('Installing the "cortix" package...')
    print('')
    !pip install cortix
    import cortix

from cortix.src.cortix_main import Cortix 
from cortix.src.network import Network

from cortix.examples.bwr.bwr import BWR
from cortix.examples.bwr.turbine import Turbine
from cortix.examples.bwr.condenser import Condenser
/home/nbuser/anaconda3_420/lib/python3.5/site-packages/matplotlib/font_manager.py:281: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  'Matplotlib is building the font cache using fc-list. '
Installing the "cortix" package...

Collecting cortix
  Downloading https://files.pythonhosted.org/packages/3e/49/94eef2ab872c1935884a1c2bf9b831e6504fd1d5a7d7aff090681b40da98/cortix-1.1.6-py3-none-any.whl (497kB)
     |████████████████████████████████| 501kB 5.2MB/s eta 0:00:01
Requirement already satisfied: scipy in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from cortix) (1.1.0)
Requirement already satisfied: numpy in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from cortix) (1.17.3)
Requirement already satisfied: matplotlib in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from cortix) (2.1.1)
Requirement already satisfied: pandas in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from cortix) (0.19.2)
Requirement already satisfied: pytest in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from cortix) (2.9.2)
Requirement already satisfied: graphviz in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from cortix) (0.13.2)
Requirement already satisfied: six>=1.10 in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from matplotlib->cortix) (1.11.0)
Requirement already satisfied: python-dateutil>=2.0 in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from matplotlib->cortix) (2.8.1)
Requirement already satisfied: pytz in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from matplotlib->cortix) (2016.6.1)
Requirement already satisfied: cycler>=0.10 in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from matplotlib->cortix) (0.10.0)
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from matplotlib->cortix) (2.1.4)
Requirement already satisfied: py>=1.4.29 in /home/nbuser/anaconda3_420/lib/python3.5/site-packages (from pytest->cortix) (1.4.31)
Installing collected packages: cortix
Successfully installed cortix-1.1.6
WARNING: You are using pip version 19.3.1; however, version 20.0.2 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-a4c304b821ca> in <module>()
     17 
     18 from cortix.examples.bwr.bwr import BWR
---> 19 from cortix.examples.bwr.turbine import Turbine
     20 from cortix.examples.bwr.condenser import Condenser

ImportError: No module named 'cortix.examples.bwr.turbine'
In [2]:
# Setup parameters
[35416] 2019-10-23 14:42:49,868 - cortix - INFO - Created Cortix object 
_____________________________________________________________________________
                             L A U N C H I N G                               
_____________________________________________________________________________
      ...                                        s       .     (TAAG Fraktur)
   xH88"`~ .x8X                                 :8      @88>
 :8888   .f"8888Hf        u.      .u    .      .88      %8P      uL   ..
:8888>  X8L  ^""`   ...ue888b   .d88B :@8c    :888ooo    .     [email protected]  @88R
X8888  X888h        888R Y888r ="8888f8888r -*8888888  [email protected]  ""Y888k/"*P
88888  !88888.      888R I888>   4888>"88"    8888    888E`    Y888L
88888   %88888      888R I888>   4888> "      8888      888E      8888
88888 `> `8888>     888R I888>   4888>        8888      888E      `888N
`8888L %  ?888   ! u8888cJ888   .d888L .+    .8888Lu=   888E   .u./"888&
 `8888  `-*""   /   "*888*P"    ^"8888*"     ^%888*     888&  d888" Y888*"
   "888.      :"      "Y"          "Y"         "Y"      R888" ` "Y   Y"
     `""***~"`                                           ""
                             https://cortix.org                              
_____________________________________________________________________________
In [3]:
# Create the application network

community = Community(n_groups=n_groups, non_offender_adult_population=100, offender_pool_size=0)
city_net.module(community)
community.end_time = end_time
community.time_step = time_step
community.show_time = (True,10*const.day)
community.save = True

prison = Prison(n_groups=n_groups, pool_size=0)
city_net.module(prison)
prison.end_time = end_time
prison.time_step = time_step
prison.save = True

parole = Parole(n_groups=n_groups, pool_size=0)
city_net.module(parole)
parole.end_time = end_time
parole.time_step = time_step
parole.save = True

adjudication = Adjudication(n_groups=n_groups, pool_size=0)
city_net.module(adjudication)
adjudication.end_time = end_time
adjudication.time_step = time_step
adjudication.save = True

jail = Jail(n_groups=n_groups, pool_size=0)
city_net.module(jail)
jail.end_time = end_time
jail.time_step = time_step
jail.save = True

arrested = Arrested(n_groups=n_groups, pool_size=0)
city_net.module(arrested)
arrested.end_time = end_time
arrested.time_step = time_step
arrested.save = True

probation = Probation(n_groups=n_groups, pool_size=0)
city_net.module(probation)
probation.end_time = end_time
probation.time_step = time_step
probation.save = True

city_net.connect( prison, parole,'bidirectional' )
city_net.connect( adjudication, prison )
city_net.connect( jail, prison )
city_net.connect( adjudication, jail )
city_net.connect( arrested, jail )
city_net.connect( arrested, adjudication )
city_net.connect( arrested, probation )
city_net.connect( probation, jail )
city_net.connect( adjudication, probation )
city.save = True

city_net.connect( prison, community )
city_net.connect( jail, community )
city_net.connect( arrested, community, 'bidirectional' )
city_net.connect( parole, community )
city_net.connect( adjudication, community )
city_net.connect( probation, community )

Verify the network connectivity

In [4]:
# View the Cortix network created

city_net.draw()
Out[4]:
network-0 0 Community 5 Arrested 0->5 1 Prison 1->0 2 Parole 1->2 2->0 2->1 3 Adjudication 3->0 3->1 4 Jail 3->4 6 Probation 3->6 4->0 4->1 5->0 5->3 5->4 5->6 6->0 6->4

Run network simulation

In [5]:
# Run the simulation!

city.run()
[35416] 2019-10-23 14:42:50,158 - cortix - INFO - Launching Module <cortix.examples.city_justice.community.Community object at 0x1514923710>
[35416] 2019-10-23 14:42:50,165 - cortix - INFO - Launching Module <cortix.examples.city_justice.prison.Prison object at 0xf146bbeb8>
[35416] 2019-10-23 14:42:50,168 - cortix - INFO - Community::time[d] = 0.0
[35416] 2019-10-23 14:42:50,172 - cortix - INFO - Launching Module <cortix.examples.city_justice.parole.Parole object at 0x1514923780>
[35416] 2019-10-23 14:42:50,182 - cortix - INFO - Launching Module <cortix.examples.city_justice.adjudication.Adjudication object at 0x15149236d8>
[35416] 2019-10-23 14:42:50,189 - cortix - INFO - Launching Module <cortix.examples.city_justice.jail.Jail object at 0x1514923748>
[35416] 2019-10-23 14:42:50,201 - cortix - INFO - Launching Module <cortix.examples.city_justice.arrested.Arrested object at 0x15148fce80>
[35416] 2019-10-23 14:42:50,212 - cortix - INFO - Launching Module <cortix.examples.city_justice.probation.Probation object at 0x1514923d68>
[35416] 2019-10-23 14:42:50,638 - cortix - INFO - Community::time[d] = 10.0
[35416] 2019-10-23 14:42:50,925 - cortix - INFO - Community::time[d] = 20.0
[35416] 2019-10-23 14:42:51,299 - cortix - INFO - Community::time[d] = 30.0
[35416] 2019-10-23 14:42:51,660 - cortix - INFO - Community::time[d] = 40.0
[35416] 2019-10-23 14:42:52,186 - cortix - INFO - run()::Elapsed wall clock time [s]: 2.32

Results inspection through Cortix

In [6]:
community     = city_net.modules[0]
prison             = city_net.modules[1]
parole             = city_net.modules[2]
adjudication   = city_net.modules[3]
jail                   = city_net.modules[4]
arrested          = city_net.modules[5]
probation        = city_net.modules[6]
quant_names  = {'Prison':'fpg','Parole':'feg','Adjudication':'fag',
                            'Jail':'fjg','Arrested':'frg','Probation':'fbg','Community':'f0g'}

total_num_params = 0
In [7]:
'''Inspect Data Function'''
def inspect_module_data(module,quant_names):
    population_phase = module.population_phase
    if isinstance(quant_names,dict):
        (fxg_quant, time_unit) = population_phase.get_quantity_history(quant_names[module.name])
    elif isinstance(quant_names,str):
        (fxg_quant, time_unit) = population_phase.get_quantity_history(quant_names)
    fxg_quant.plot( x_scaling=1/const.day, x_label='Time [day]', y_label=fxg_quant.name+' ['+fxg_quant.unit+' %]')
    plt.grid()
    plt.show()

    # number of parameters in the prison model
    n_params = (len(population_phase.GetActors())-1)*n_groups
    print('# parameters = ',n_params)
    return n_params

Prison module

In [8]:
total_num_params += inspect_module_data(prison,quant_names)
# parameters =  1800

Parole module

In [9]:
total_num_params += inspect_module_data(parole,quant_names)
# parameters =  1800

Adjudication module

In [10]:
total_num_params += inspect_module_data(adjudication,quant_names)
# parameters =  3600

Jail module

In [11]:
total_num_params += inspect_module_data(jail,quant_names)
# parameters =  1800

Arrested module

In [12]:
total_num_params += inspect_module_data(arrested,quant_names)
# parameters =  3600

Probation module

In [13]:
total_num_params += inspect_module_data(probation,quant_names)
# parameters =  1800

Community module

Offenders removed (negative sign) from the general non-offender community.

In [14]:
total_num_params += inspect_module_data(community,quant_names)
# parameters =  2250
In [15]:
total_num_params += inspect_module_data(community,'f0g_free')
# parameters =  2250
In [16]:
'''Total number of unknowns and parameters'''

print('total number of unknowns   =', n_groups*len(city_net.modules))
print('total number of parameters =', total_num_params)
total number of unknowns   = 3150
total number of parameters = 18900