• Không có kết quả nào được tìm thấy

Modeling and Simulation in Python

N/A
N/A
Nguyễn Gia Hào

Academic year: 2023

Chia sẻ "Modeling and Simulation in Python"

Copied!
247
0
0

Loading.... (view fulltext now)

Văn bản

Permission is granted to copy, distribute, transmit and adapt this work under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License: http://modsimpy.com/license. The source and LATEX code for this book are available from https://github.com/AllenDowney/ModSimPy.

Can modeling be taught?

How much programming do I need?

How much math and science do I need?

Getting started

Installing Python

Copying my files

If you download the Zip file, you can skip the rest of this section, which explains how to use Git. Before you download these files, I suggest you copy my repository on GitHub, which is called forking.

Running Jupyter

The goal of the modeling process is to find the model that best suits its purpose (prediction, explanation, or design). As an example, suppose someone asks you why the Earth's orbit is elliptical.

Computation

In this chapter's workbook, you will see another model, which assumes that the acceleration is constant until the penny reaches terminal velocity. Next, we calculate the time it takes for the penny to fall 381 m, the height of the Empire State Building.

Defining functions

When you call the function, it executes the statements in the body, which update the variables of the bikeshare object; you can check it by displaying the new state. Another benefit is that the name you give the function documents what it does, which makes programs more readable.

Print statements

If statements

Parameters

The advantage of using parameters is that you can call the same function many times, providing different arguments each time. Adding parameters to a function is called generalization because it makes the function more general, that is, less specialized.

For loops

TimeSeries

When the loop exits, the results contain 10 timestamps, from 0 to 9, and the number of bikes at Olin at the end of each time step. TimeSeries is a specialized version of Series, which is defined by Pandas, one of the libraries we'll be using a lot.

Plotting

Before moving on, take a minute to review the model from the previous chapters. The model does not check if a bike is available, so it is possible for the number of bikes to be negative (as you may have noticed in some of your simulations).

Figure 2.1: Simulation of a bikeshare system showing number of bikes at Olin over time.
Figure 2.1: Simulation of a bikeshare system showing number of bikes at Olin over time.

More than one State object

As long as there is only one state facility, that's fine, but what if there are more than one bike sharing system in the world. So we can simulate different bike sharing systems, or run multiple simulations of the same system.

Documentation

Negative bikes

Comparison operators

Metrics

We've seen several functions that return values; for example, when you run sqrt, it returns a number that you can assign to a variable. For example, if you run step, it updates a State object, but returns no value.

Two kinds of parameters

Loops and arrays

Each time the loop is run through, the next number in the sequence is assigned to the loop variable, i. The arguments specify where the sequence should start and stop and how many elements it should contain.

Sweeping parameters

Incremental development

In this chapter and the next, we use the tools from the previous chapters to explain world population growth since 1950 and create predictions for the next 50–. For background on world population growth, check out this video from the American Museum of Natural History http://modsimpy.com/human.

Plotting

Constant growth model

In the next few sections, I demonstrate this process, starting with simple models and gradually improving them. Although there is some skew in the plotted estimates, world population growth since 1960 appears to have been nearly linear.

Simulation

But first, we can improve the code from the previous chapter by encapsulating it in a function and using System objects. State objects contain state variables, which represent the state of the system, which are updated during a simulation.

Figure 5.2: Estimates of world population, 1950–2016, and a constant growth model.
Figure 5.2: Estimates of world population, 1950–2016, and a constant growth model.

Proportional growth model

Now we can choose the values ​​of birth_rate and death_rate that best fit the data. The proportional model fits the data well from 1950 to 1965, but not so well thereafter.

Factoring out the update function

Rather than repeating identical code, we can separate the things that change from the things that don't. When we call run_simulation, the second parameter is a function, such as update_func1, that calculates the population for the next year.

Figure 6.1: Estimates of world population, 1950–2016, and a proportional model.
Figure 6.1: Estimates of world population, 1950–2016, and a proportional model.

Combining birth and death

When the population is less than 3-4 billion, net growth is proportional to the population, as in the proportional model. Just below 14 billion there is a point where net growth is 0, meaning the population does not change.

Figure 7.1: Estimates of world population, 1950–2016, and a quadratic model.
Figure 7.1: Estimates of world population, 1950–2016, and a quadratic model.

Equilibrium

At this point, the birth and death rates are equal, so the population is unequal.

Dysfunctions

The use of "projection" leaves open the possibility that there are important things in the real world that are not captured in the model. It also suggests that even if the model is good, the parameters we estimate based on the past may be different in the future.

Figure 8.1 shows the result, with a projection until 2250. According to this model, population growth will continue almost linearly for the next 50–100 years, then slow over the following 100 years, approaching 13.9 billion by 2250.
Figure 8.1 shows the result, with a projection until 2250. According to this model, population growth will continue almost linearly for the next 50–100 years, then slow over the following 100 years, approaching 13.9 billion by 2250.

Comparing projections

So if we want to know x100 and don't care about the other values, we can calculate it with a multiplication and an addition. There is no analytical solution to this equation, but we can approximate it to a differential equation and solve it, which we will do in the next section.

Figure 8.2: Projections of world population generated by the U.S. Census Bureau, the United Nations, and our quadratic model.
Figure 8.2: Projections of world population generated by the U.S. Census Bureau, the United Nations, and our quadratic model.

Differential equations

Since K is an arbitrary constant, exp(K) is also an arbitrary constant, so we can write . If you want to see this derivation in more detail, you might like this video: http://modsimpy.com/khan1.

Analysis and simulation

With analysis we can sometimes accurately and efficiently calculate a value that we could approach less efficiently with simulation alone. For example, in Figure 7-2, we can see net growth going to zero near 14 billion, and we could estimate carrying capacity using a numerical search algorithm (more on that later).

Analysis with WolframAlpha

Analysis with SymPy

Python does not attempt to perform numerical addition; rather, it creates a new Symbol that represents the sum of t and 1. If we now write f(t), we get an object that represents the evaluation of a function, f, at a value, t.

Differential equations in SymPy

Solving the quadratic growth model

Summary

Instead of returning the final state of the system, it returns a Time Series containing the state as it changes over time. The update function takes the current state of the system, time and system object as parameters and returns the new state.

Under the hood

Each DataFrame contains three objects: index is an array of labels for the rows, column is an array of labels for the columns, and values ​​is a NumPy array containing the data. The type table2.valuesis ndarray, which is the primary array type provided by NumPy; the name indicates that the array is "n-dimensional"; that is, it can have an arbitrary number of dimensions.

Figure 10.1: The elements that make up a DataFrame and a Series.
Figure 10.1: The elements that make up a DataFrame and a Series.

One queue or two?

As with the bike sharing model, we assume that a customer is equally likely to arrive at each time step. The value of λ probably varies from day to day, so we will have to consider a number of possibilities.

Figure 10.2: One queue, one server (left), one queue, two servers (middle), two queues, two servers (right).
Figure 10.2: One queue, one server (left), one queue, two servers (middle), two queues, two servers (right).

Predicting salmon populations

In the archive for this book you will find a notebook called queue.ipynb that contains some code to get started and instructions.

Tree growth

Depending on the quality of the site, the trees can grow faster or slower. So each curve is identified by a "site index" that indicates the site's quality.

Figure 10.3: Site index curves for tree growth.
Figure 10.3: Site index curves for tree growth.

The Freshman Plague

In this chapter, we develop a model of an epidemic as it spreads in a susceptible population and use it to evaluate the effectiveness of possible interventions.

The SIR model

The SIR equations

Sub-fund models are often represented visually with stock and flow charts (see http://modsimpy.com/stock). The widget in the center of the arrows represents a valve that controls the flow rate; the diagram shows the parameters controlling the valves.

Implementation

The update function

The assignment does exactly what we want: it assigns the three values ​​of the State object to the three variables, in order. You may notice that this version of update_func does not use any of its parameters,t.

Running the simulation

Collecting the results

The proportion of the infected population is never very high, but it adds up.

Figure 11.2: Time series for S, I, and R over the course of 98 days.
Figure 11.2: Time series for S, I, and R over the course of 98 days.

Now with a TimeFrame

In the SIR model, we may want to know the time to the peak of the outbreak, the number of people who are sick at the peak, the number of students who will still be sick at the end of the semester, or the total number of students who become sick at any given time. As an example, I will focus on the latter: the total number of sick students.

Immunization

If the series is a time series, the label you get from idxmax is a time or date. The result is a SweepSeries object that maps from each immunization rate to the resulting proportion of students ever infected.

Figure 12.1: Time series for S, with and without immunization.
Figure 12.1: Time series for S, with and without immunization.

Hand washing

How should we quantify the impact of the money we spend on a handwashing campaign? Now we need to model the relationship between the money we spend and the effectiveness of the campaign.

Figure 12.3: Fraction of the population infected as a function of hand-washing campaign spending.
Figure 12.3: Fraction of the population infected as a function of hand-washing campaign spending.

Optimization

For each number of doses, we calculate the fraction of students we can immunize, fraction and the remaining budget we can spend on the campaign. As we increase the number of doses, we have to cut campaign spending, which seems to make things worse.

Figure 12.4 shows the result. If we buy no doses of vaccine and spend the entire budget on the campaign, the fraction infected is around 19%
Figure 12.4 shows the result. If we buy no doses of vaccine and spend the entire budget on the campaign, the fraction infected is around 19%

Sweeping gamma

When gamma is low, the recovery rate is low, meaning people are contagious for longer.

SweepFrame

This figure suggests that there may be a correlation between beta and gamma that determines the outcome of the model. The figures in Section 13.3 suggest that there is a relationship between the SIR model parameters, beta and gamma, which determine the outcome of the simulation, the proportion of students who are infected.

Figure 13.3: Fraction of students infected as a function of the parameter gamma, for several values of beta.
Figure 13.3: Fraction of students infected as a function of the parameter gamma, for several values of beta.

Exploring the results

Describing physical systems using dimensionless parameters is often a useful move in the modeling and simulation game. In the notebook for this chapter, you can explore the first option as an exercise.

Contact number

By convention it is usually labeled R0, but in the context of the SIR model this notation is confusing, so we will use c instead. The results in the previous section show that there is a relationship between the total number of infections.

Analysis and simulation

Recall from Section 10.2 that an Array object contains an index and a corresponding sequence of values. When the contact number is small, the early progress of the epidemic depends on details of the scenario.

Figure 14.2: Total fraction infected as a function of contact number, showing results from simulation and analysis.
Figure 14.2: Total fraction infected as a function of contact number, showing results from simulation and analysis.

Estimating contact number

So far the systems we have studied have been physical in the sense that they exist in the world, but they have not been physics in the sense that physics classes are usually about. In the next few chapters we're going to do some physics, starting with thermal systems, that is, systems where the temperature of objects changes as heat is transferred from one to another.

The coffee cooling problem

To use this data and answer the question, we need to know something about temperature and heat, and we need to make some modeling decisions.

Temperature and heat

Heat transfer

Fluid flow can be caused by external action, such as stirring, or by internal differences in temperature. Convection can be a complex subject, as it often depends on details of fluid flow in three dimensions.

Newton’s law of cooling

Convection: When particles in a gas or liquid flow from place to place, they carry heat energy with them. Radiation: When particles in an object move due to thermal energy, they emit electromagnetic radiation.

Implementation

We can adjust r and find the right value by trial and error, but we will see a better way in the next chapter. In general, we do not know the value of r, but we can use measurements to estimate it.

Mixing liquids

Assuming there are no chemical reactions that either produce or consume heat, the total thermal energy of the system is the same before and after mixing; in other words, thermal energy is conserved. If the temperature of the first fluid is T1, the temperature of the second fluid is T2, and the final temperature of the mixture is T, the heat transfer to the first fluid is C1(T − T1) and the heat transfer to the second fluid is C2(T − T2), where C1 and C2 are the thermal masses of the fluids.

Mix first or last?

Optimization

The end temperature is maximized whent_add=30, so adding the milk at the end is optimal. In that case, it's better to add the milk at the coffee shop, or wait until I get to the office.

Figure 16.2 shows the result. Again, note that this is a parameter sweep, not a time series
Figure 16.2 shows the result. Again, note that this is a parameter sweep, not a time series

Analysis

In this chapter, we implement one of the most widely used pharmacokinetic models: the so-called minimal model of glucose and insulin in the bloodstream. Prolonged, severe hyperglycemia is. (see http://modsimpy.com/cdc) to define the symptoms of diabetes, a serious disease that affects nearly 10% of the US population.

The glucose minimal model

X is the concentration of insulin in the tissue fluid as a function of time, and dX/d is its rate of change. I is the blood insulin concentration versus time, which is taken as input to the measurement-based model.

Data

Interpolation

Next, we will see how to find the parameters that produce the series that best fits the data. The graph above shows the simulated glucose levels from the model along with the measured data.

Figure 18.1 shows the results. The top plot shows simulated glucose levels from the model along with the measured data
Figure 18.1 shows the results. The top plot shows simulated glucose levels from the model along with the measured data

Solving differential equations

The ModSimSeries from run_ode_solver that we assign to details contains information about how the solver is running, including a success code and diagnostic message. In Chapter 16, we used root_bisect to find the value of a parameter that produces a particular result.

Figure 18.2: Results from Euler’s method (simulation) and Ralston’s method (ODE solver).
Figure 18.2: Results from Euler’s method (simulation) and Ralston’s method (ODE solver).

The insulin minimal model

Low-Pass Filter

A filter is "low-pass" if it allows low-frequency signals to pass unchanged from Vin to Vout, but it reduces the amplitude of high-frequency signals. By applying the laws of circuit analysis, we can derive a differential equation that describes the behavior of this system.

Thermal behavior of a wall

Newton's law may not seem like a differential equation until we realize that acceleration, a, is the second derivative of position, y, with respect to time, t. If mass depends on time, position, or velocity, we must use a more general form of Newton's law (see http://modsimpy.com/varmass).

Dropping pennies

It is not used in this slope function because none of the factors of the model are time dependent (see section 18.1). The rest of the function is a simple translation of the differential equations, with the substitution a=−g, indicating that the acceleration due to gravity is in the direction of decreasing y.

Figure 20.1: Height of the penny versus time, with no air resistance.
Figure 20.1: Height of the penny versus time, with no air resistance.

Events

In this context, the reference area is the projected frontal area, that is, the visible area of ​​the object as seen from a point on its line of travel (and far away). Cd is the drag coefficient, a dimensionless quantity that depends on the shape of the object (including length, but not frontal area), its surface properties and how it interacts with the fluid.

Implementation

In the next chapter, we will use vector objects to keep track of the direction of the forces and add them together in a less error-prone way. In the previous chapter, we modeled objects moving in one dimension, with and without drag.

Vectors

In the Cartesian plane, the angle 0 rad is east, and the angle πrad is west. It can be difficult and error-prone to use both units in the same program.

Simulating baseball flight

Another way to represent the direction of A is a unit vector, which is a vector with magnitude 1 pointing in the same direction as A. We can do the same thing using the hat function, so named because unit vectors are conventionally decorated with a hat, like this: ˆA.

Trajectories

Finding the range

Finishing off the problem

Implementation

Analysis

Moment of inertia

Teapots and turntables

Estimating friction

Bungee jumping

Bungee dunk revisited

Spider-Man

Kittens

Simulating a yo-yo

Hình ảnh

Figure 2.1: Simulation of a bikeshare system showing number of bikes at Olin over time.
Figure 5.2: Estimates of world population, 1950–2016, and a constant growth model.
Figure 6.1: Estimates of world population, 1950–2016, and a proportional model.
Figure 7.1: Estimates of world population, 1950–2016, and a quadratic model.
+7

Tài liệu tham khảo

Tài liệu liên quan

You are now ready to implement the circuit. The steps in creating the circuit will be as follows. Power will be provided to the breadboard. A switch will be inserted into