Dynamic systems -------------------------------------------- - - PDF document
Lesson 2 University of Bergamo Engineering and Management for Health FOR CHRONIC DISEASES MEDICAL SUPPORT SYSTEMS LESSON 2 Overview of the main methods we will adopt during the course: differential methods; stochastic methods; optimization
Ettore Lanzarone March 11, 2019
ORDINARY DIFFERENTIAL EQUATION Let us consider a function y(t). A ordinary differential equation (ODE) is: The goal is to find the solution y(t) that respect the equation tI. The solution is also called “integral” of the equation (as it is usually computed through integrals). The highest order of derivative included in the equation is the order of the differential equation.
PARTIAL DIFFERENTIAL EQUATION Let us consider a function y(t,x1,…,xm). A partial differential equation (PDE) is: (limited to the second order) are the partial derivatives with respect to one of the independent variables (computed fixing the values of the other variables and considering them as constant values)
are the second partial derivatives with respect to one of the independent variables Second partial derivatives can computed with respect to two variables (first time derivative with respect to one variable, second time derivative with respect to the other variable): MIXED PARTIAL DERIVATIVES. Schwarz's theorem: mixed partial derivatives are the same while changing the order of the derivatives, e.g.;
Let us consider a simple ODE: where f(t) is a known function. The solution is: There are infinite solutions that depend on the value of coefficient k. General validity: an ODE of order n has infinite solutions that depend on n coefficients. The set these solutions is named GENERAL INTEGRAL of the ODE.
In practical application, we are interested in choosing one of the solutions. How? Usually, initial condition:
The single solution that respects the condition is named PARTICULAR INTEGRAL. In general, for an ODE of order n we need n initial conditions: y(t0), y’(t0), y’’(t0), …, y(n-1)(t0)
Other alternatives are also possible, e.g., to impose that the solution passes through several points. As for PDEs, the idea is similar. We have to include BOUNDARY CONDITIONS. See the example for y(t,x)
Possible approaches to solve a differential equation:
Most powerful approach: it allows to get the analytical expression of the solution. Limitation: analytical solutions are available only for a limited number of equation classes
The idea is to draw a qualitative trend of the solution using the information that can be qualitatively extracted from the solution (e.g., the signum of the derivative to get local minima and maxima)
Flexible approach: it can be always applied. Limitations: computational time to get the solution, solution given in terms of points without getting the analytical expression; convergence of the numerical approach to be verified
EXAMPLE OF ANALYTICAL SOLUTION SEPARABLE VARIABLES 1. Constant solutions:
EXAMPLE OF ANALYTICAL SOLUTION SEPARABLE VARIABLES 2. Other solutions:
EXAMPLE OF ANALYTICAL SOLUTION SEPARABLE VARIABLES 2. Other solutions:
EXAMPLE OF ANALYTICAL SOLUTION SEPARABLE VARIABLES Constant solutions: Other solutions:
EXAMPLE OF ANALYTICAL SOLUTION LINEAR FIRST ORDER EQUATION a(t) and f(t) are known continuous functions. The solution is:
These are just examples to show analytical solutions. Probably none of the problems in medicine and health care can be solved in this way. We will focus on the numerical approaches.
Possible approaches to solve a differential equation:
Most powerful approach: it allows to get the analytical expression of the solution. Limitation: analytical solutions are available only for a limited number of equation classes
The idea is to draw a qualitative trend of the solution using the information that can be qualitatively extracted from the solution (e.g., the signum of the derivative to get local minima and maxima)
Flexible approach: it can be always applied. Limitations: computational time to get the solution, solution given in terms of points without getting the analytical expression; convergence of the numerical approach to be verified
Numerical methods are based on the discretization of the continuous independent variable x or t in y(x) or y(t), respectively. Instead of considering the entire function y(x) we consider samples at values xi.
We consider a first order ODE: We get estimates of the solution at different base points:
0)
The solution at the next point is written in terms of the Taylor polynomial: Polynomial centered in x0 and evaluated in x=x0+h
n n x x n n x x x x n k n k x x k k
2 2 2
The simpler approximation is to cut at the first order: We may express the derivative based on the ODE: Thus:
x x
0
x x
The approach is iterative: where y(x0+h) is known from the previous step where y(x0+2h) is known from the previous step
In a compact way: with:
1
i i i i i
GEOMETRICAL INTERPETATION
x0 x y0
INITIAL POINT (x0,y0)
GEOMETRICAL INTERPETATION
h x0 x1 x y0
y1=y0+f(x0,y0)h
Slope=f(x0,y0)
f(x0,y0)h
y1
GEOMETRICAL INTERPETATION
x0 x1 x2 x y0
y1=y0+f(x0,y0)h
y1 h
y2=y1+f(x1,y1)h
New slope=f(x1,y1) y2
hf(x1,y1)
Example with h=0.01
2
2 2 2 2 3 2 1 1 1 2 2 1 1
i i i i
Example True value is known because for this example also the analytical solution exists (separable variables)
Example The solution is given by the coordinates of the discrete points and can be easily plotted.
The problem is to determine an integration step h that is small enough to get convergence. Suggestion:
In some cases where the slopes of the function y(x) are high the step h can be really small or the convergence is never guaranteed.
In some cases where the slopes of the function y(x) are high the step h can be really small
Use something else with respect to the Euler method: 1. Higher order methods: do not limit to the first order in the Taylor polynomial that approximates f(x,t) 2. Add additional points between xi and xi+1=xi+h RUNGE-KUTTA METHODS
Runge-Kutta denotes a family of method. Among them, the RK4 is widely adopted. with Instead of using only one slope k=k1, a set of four slopes are considered to better approximate the solution.
Graphical representation of the slopes The considered slope is then a wighted average of them.
Numerical approaches can be extended to study systems of differential equations:
An ODE with order n larger than one can be rewritten as a system of n ODE equations with order 1. Example for second order: We define two new variables: So the equation becomes:
MATLAB
SIMNON Software tool to easily solve ODEs and systems of ODEs. You can find a free version within the course material. Let’s use the program! 1. Write the system in a text file *.t
SIMNON Software tool to easily solve ODEs and systems of ODEs. You can find a free version within the course material. Let’s use the program! 1. Write the system in a text file *.t 1. Run the model using the command line
EXPORT store < store Results are stored in a file store.d This command generates a text file with the discrete points of the solution, named store.t
Note:
The integration step h can be optimized by the solver in
computational performance while keeping the convergence of the solution.
It is possible to decide the parameters of the simulation. In SIMNON, instead of the simple command line “SIMU 0 ***” it is possible to configure the parameters with a dedicated window. CHOOSE AUTO OR A FIXED STEP
It is also possible to choose the integration algorithm
Optimization means to find the maximum and minimum of a function over a given set (entire domain of the function or a subdomain). If the set coincides with the domain it is UNCONSTRAINED OPTIMIZATION;
Function y(x) - unconstrained
Maximum and minimum points are searched among those with y’(x)=0
For example: y=|x| has a minimum in x=0, which is not determined with y’(x)=0
Function y(x) - constrained
The maximum is on the border of the interval if the set is close [a,b] If the set is open (a,b) there is no maximum but only an upper extreme
Function y(x) - constrained
Easy to see but important. It will be at the basis of the simplex method for linear programming models. Maximum and minimum on the border of the interval if the set is close [a,b] If the set is open (a,b) only upper and lower extremes
Function z=f(x,y) – unconstrained
Maximum and minimum points are searched among those with null gradient.
in these points:
Hessian is a symmetric matrix (Schwarz's theorem: mixed partial derivatives are equal)
Hessian cases: 1. Positive determinant and positive first element MINIMUM 2. Positive determinant and negative first element MAXIMUM 3. Negative determinant SADDLE POINT 4. Null determinant NO INFORMATION Similar to y’’(0)=0 for functions of 1 variable
SADDLE POINT
Function z=f(x,y) - constrained
The border is a close line expressed as a function g(x,y)=0 Two possible approaches: 1. Parametrization of the border 2. Lagrange multipliers
Parametrization of the border The border g(x,y)=0 can be rewritten as: The problem is to find maximum and minimum for a function of only one variable.
Lagrange multipliers method The following function is built based on z=f(x,t) and g(x,y)=0 The candidate maxima and minima are those with a null gradient of (x,y,)
Then, computing the value of z in these points, maxima and minima are identified.
In case the function is linear and the border is piecewise composed by linear equations (see for example the figure), the minimum and the maximum of the function in the domain are
In this case, the problem can be rewritten as a so-called linear programming model. Or alternatively with maximize All inequality are with “=“ to include the border (maximum and minimum are on the border)
This problem can be solved via a set of algorithms, e.g., the SIMPLEX ALGORITHM. Commercial solvers for this type of problems. For example, CPLEX couples simplex and some branch and bound / branch and cut algorithms.
Free version for students and researchers available: www.ibm.com/it-it/products/ilog-cplex-optimization-studio
For problems with long computational times, CPLEX provides intermediate solutions associated with an optimality gap. When gap is 0, the optimal solution is found. It is possible to stop the execution after a given computational time or when a solution associated with a given gap is obtained. What is the gap? Estimation of the difference between the objective function of the current solution and that of the best possible one.
GAP = estimation of the difference between the objective function of the current solution and that of the best possible one. For example, in case of problems with integer variables, a relaxed problem is solved using continuous variables with respect to integer ones. Then the gap is based on this relaxed solution. Obviously that solution is not admissible for the original problem. The gap is updated because: