Modeling Basics
Compartment Models
Dimensional Analysis
Stochastic Modeling
ME 635: Modeling and Simulation
Kishore Pochiraju
Model
Physics Math
Simulation
Computation
Numerical Methods
Basic Concepts: Differentiability
• A function is differentiable if there is a matrix [A] such that:
Linearization:
Jacobian
Solution to Linear Simultaneous Equations
{x} is non-trivial if:
Solving by Cramer’s Rule
Solving by Gaussian Elimination
The general procedure for Gaussian Elimination can be
summarized in the following steps:
– Write the augmented matrix for
the system of linear equations.
– Use elementary row operations
on the augmented matrix [A|b]
to transform A into upper
triangular form.
• If a zero is located on the diagonal,
switch the rows until a nonzero is in
that place.
• If you are unable to do so, stop; the
system has either infinite or no
solutions.
– Use back substitution to find
the solution of the problem.
Computational Performance
• Compare number of required operations for solving linear
system equations.
Time for Operations
Memory Issues
• Single Precision Float data variable requires 4 Bytes (32
Bits):
IEEE 754 standard binary32
• Double Precision 8 Bytes – 64 bits
Storing the Matrix K.
[K] has N2 data values
– Full Storage: 8*N2 Bytes (Double precision)
– Symmetric: 8 * N(N+1)/2 Bytes
– Sparse & Banded – Less Storage
Sparse Matrix (non-zeros are dark) Banded Matrix
Indirect / Iterative Methods reduce Memory Burden
Why are eigenvalues important?
Compartment Models
Input Rate of change in Output
rate compartment rate
Water level
Water outflow
Water inflow
Reservoir
Typical Compartment Model
Room
Heat in Heat out
Temperature (q)
Q_in Q_out
Rate of Change in Compartment = input rate – output rate
Typical Solutions to Population Compartments
• Sigmoid functions -
Appropriate for
population growth,
diffusion, etc.
An Initial Value Problem
Rate of change Linear ODE
with time
Solution:
What is the initial value?
What is the equilibrium value?
Time constants of ODES
• Consider:
• At steady state : x=x0 : f(x0) = 0
• Near equilibrium, if f is differentiable:
Time constants are related to eigenvalues of A
Eigen expansion of A if
diagonalization is possible
What are stiff ODEs
• ODE is stiff around the equilibrium solution if its time
constants that range over several orders of magnitude.
• Stiffness is important when choosing/designing numerical
methods for solving ODES
Linear Regression
Steps:
• Assume a function with known functions and unknown
coefficients
• Solve the unknown coefficients by minimization.
Model from Observations
• Objective: We observed 10 states of a system, i.e. if y = f(x),
we have 10 values of (x,y),
Use a model to estimate the value of y for a value of x
where there is no observation
X Y
1 0.3415069
3 1.0920367
4 1.4187925
7 1.6292369
8 1.3082736
8.5 1.0623416
9 0.7602542
9.1 0.693305
10 0.0002519
Observations/Samples & Interpolation
Observations
Piece-wise linear
Cubic
linear
• Exact Model
Interpolation/Extrapolation
• Interpolation:
data to be found
are within the
Exact
range of observed
data
• Extrapolation:
data to be found
are beyond the
range of
Cubic observation
data. (not reliable*)
Extrapolation should
With unknown/incorrect functional be avoided whenever
forms extrapolation rarely works! possible
The Regression or Curve Fit Problem
• Definition:
– Given N observations (xi, yi) with a curve/surface with M
unknown coefficients, aj, with N > M.
• Strategy:
– Find the coefficients a, such the to square error of distance
between the model and observations is minimized.
Numerical Example
Error: ei
How to minimize?
Say, the model consists of M linear coefficients and known functions:
Now we have N observations. We can write: (error will be a function of the
unknown coefficients)
Coefficients that minimize the error are:
Example: If
Y = a0 + a1 x = a0 (1) + a1 (x)
The two known functions are: (1) and (x)
Say we have several (M) observations: we
get
Best Fit {a}
• {a} that minimizes for this equation is:
Solving
Can be used for determination of constant coefficients & Any known
functions
y1 1 x1 x12 x13 1
a
1 ... b
1 ...
1 c
...
d
3
yn 1 x n xn2 xn n
32
Implementation in MATLAB
Consider the following data set:
x y
1 4.256430363
3 10.76433253
10 33.36067977
20 64.51141009 This data was created with:
24.5 77.83651186
30.9 95.81242635 Y = 1 + 2*x + 40* sin(px/100)
37.3 112.4582736
43.7 127.6191049
50.1 141.1998026
56.5 153.1689124
62.9 163.5599109
69.3 172.4698563
75.7 180.0554788
82.1 186.5268648
88.5 192.1389938
100 201
Read the data file into Matlab
• [x,y] = textread('/Users/kishore/Documents/Courses/Modeling and
SImulation -2011/Fits_Dataset-1.csv','%f,%f’)
– Textread – Reads the text files
– File format is ‘%f,%f’ = float, float , e.g: 1.0,1.1
– Reads the file into two arrays, x,y
• Lets say we want to do the best fit line into it:
• Y=Mx +C
Then
[X] = [ 1 Xi]
Construction of [X]
• temp1 = ones (length(x),1)
– Constructs a column of ones.
• X = [temp1 x]
Constructs a [length(x) x 2 ] matrix
Y = y – Same as y – output observations
The Regression or Curve Fit Problem
• Definition:
– Given N observations (xi, yi)
with a curve/surface with M
unknown coefficients, aj, with
N > M.
• Strategy:
– Find the coefficients “a”, such
the square error of distance
between the model and
observations is minimized.
Polynomial Regression
Best Fit Approximations with Least Square Methods
• What if the number of observations are greater than the
parameters in the model.
• For example: Many observations nearly linear behavior
Numerical Example
Error: ei
Implementation in MATLAB
Consider the following data set:
x y
1 4.256430363 This data was created with:
3 10.76433253
10 33.36067977 Y = 1 + 2*x + 40* sin(px/100)
20 64.51141009
24.5 77.83651186
30.9 95.81242635
37.3 112.4582736
43.7 127.6191049
50.1 141.1998026
56.5 153.1689124
62.9 163.5599109
69.3 172.4698563
75.7 180.0554788
82.1 186.5268648
88.5 192.1389938
100 201
Read the data file into Matlab
• [x,y] = textread('/Users/kishore/Documents/Courses/Modeling and
SImulation -2011/Fits_Dataset-1.csv','%f,%f’)
– Textread – Reads the text files
– File format is ‘%f,%f’ = float, float , e.g: 1.0,1.1
– Reads the file into two arrays, x,y
• Lets say we want to do the best fit line into it:
• Y=Mx +C
Then
[X] = [ 1 Xi]
Construction of [X]
• temp1 = ones (length(x),1)
– Constructs a column of ones.
• X = [temp1 x]
Constructs a [length(x) x 2 ] matrix
Y = y – Same as y – output observations
Matlab
[x,y] = textread('/Users/kishore/Documents/Courses/Modeling and SImulation
-2011/Fits_Dataset-1.csv','%f,%f')
plot(x,y,'*')
hold on
% Linear Interpolation y = mx + c
temp1 = ones (length(x),1)
X = [temp1 x];
coeffs = inv(X'*X)*(X'*y)
% Generate the line
p = [1:100]';
q = coeffs(1) + coeffs(2)*p;
plot(p,q)
Figure
Fits
Linear Quadratic
cubic Exact
Matlab Example
x = [0.9 1.5 3 4 6 8 9.5];
y = [0.9 1.5 2.5 5.1 4.5 4.9 6.3] ;
plot(x, y, 'b*-', 'LineWidth', 2, 'MarkerSize', 15);
coeff = polyfit(x, y, 1); % coeff returns 2 coefficients fitting y = a_1
* x + a_2% Get fitted values
fittedX = 0.9:0.1: 9.5;
fittedY = polyval(coeff, fittedX);
% Plot the fitted line
hold on;
plot(fittedX, fittedY, 'r-', 'LineWidth', 2);
grid
xlabel('X')
ylabel('Y')
legend('Given data' , 'Fitted data')
Non-Deterministic Fits
(Random Sample Consensus: RANSAC)
• RANSAC is an iterative method for estimating a
mathematical model from a data set that contains outliers.
• RANSAC is accomplished with the following steps
Randomly selecting a subset of the
data set.
Fitting a model to the selected subset.
Determining the number of outliers.
Repeating steps 1-3 for a prescribed
number of iterations.
Data points shown in blue, with the line of
form y = mx+c estimated using RANSAC
indicated in red.
Interpolation
• Concentrate first on the case where we believe there is no error in the data (and
round-off is assumed to be negligible).
• So we have yi=f(xi) at n+1 points x0,x1…xi,…xn: xj > xj-1
• Conceptually, interpolation consists of two stages:
– Develop a simple function g(x) that
• Approximates f(x)
• Passes through all the points xi
– Evaluate f(xt) where x0 < xt < xn
• The selection of the simple functions g(x)
– Polynomial interpolation
– Spline (polynomial) interpolation
– Least-squares (polynomial) approximation
– Trigonometric functions (periodic)
• Fourier Transform and representations
Polynomial Interpolation
• We assume the existence of a set of data points
(xi, yi), i = 0, 1, …, n obtained from a function
f(x) so that yi = f(xi), i = 0, 1, …, n.
• In general, given n+1 points, there is a unique
polynomial gn(x) of order n:
g n ( x) a0 a1 x a2 x an x
2 n
• That passes through all n+1 points
• Polynomial Interpolation approaches:
– Lagrange interpolating polynomials
– Newton’s divided difference interpolating polynomials
Lagrange interpolating polynomials
• A suitable function for interpolation f(x) is expressible as:
n
f n ( x) Li ( x) f ( xi )
i 0
The functions Li(x), i = 0, 1, …, n are chosen to satisfy
Li ( x)
n
x xk
k 0, k i xi xk
1 i j
Li ( x j ) ij
0 i j
Where the Li are weighting coefficients that are functions of x.
Linear interpolation
• Given points A and B and x2, compute y2, yielding point C
1st order Lagrange interpolating polynomial
• 1st Lagrange interpolating
polynomial may be obtained
from a weighted combination of
two linear interpolations, as
shown.
• The resulting formula based on
known points x1 and x2 and the
values of the dependent
function at those points is:
f1 ( x) L1 f x1 L2 f x2
x x2 x x1
L1 , L2
x1 x2 x2 x1
x x2 x x1
f1 ( x) f x1 f x2
x1 x2 x2 x1
Quadratic Interpolation
• One Dimensional: (-1 < r < 1), Observations at r = -1,
0, +1
• Two Dimensional: 9 Observations @ (r,s) = (-1, -1)
(0,-1), (1,-1), (0,-1) … (0,0)
s
7
4 3
6
8 9
r
1 2
5
2nd order Lagrange interpolating polynomial
The figure shows a second order case:
f 2 ( x)
x x1 x x2 f ( x0 )
x0 x1 x0 x 2
x x0 x x2
f ( x1 )
x1 x0 x1 x 2
x x0 x x1
f ( x2 )
x2 x0 x2 x1
Each of the three terms passes through
one of the data points and zero at the
other two. The summation of the three
terms must, therefore, be unique
second order polynomial f2(x) that
passes exactly through three points.
Bilinear Interpolation in Grid Square
• What about in 2D? Q11 = (x1, y1), Q12 = (x1, y2),
– Interpolate in x, then in y Q21 = (x2, y1), and Q22 =
(x2, y2).
• Example
– We know the red values
– Linear interpolation in x between
red values gives us the blue values
– Linear interpolation in y between
the blue values gives us the
answer
http://en.wikipedia.org/wiki/
Bilinear_interpolation
2D Lagrange Interpolators
Domain:
Linear functions: (1-r)(1-s)(1+r)(1+s)
r
Value at M located at coordinate (rm, sm) with values Vi known at nodes
Piece-Wise Models and Continuity
• When using Piece-wise functions, the models may or may
not be continuous
• Enforce continuity by forcing the coefficients.
@ x=0.5, both
models produce
the same value –
No discontinuity
In y – C0
continuous
@ x=0.5, the y=20x y=-16x+18
slopes (dy/dx)
are not
continues, -- NO
C1 continuity
Examples of Smoothness/Continuity
Continuity of Piece-Wise Functions
• For y=f(x) and y=g(x), C0 exists at the interface between
two models (x=x0), if the value of the y (prediction) is the
same ; f(x0) = g(x0)
• E.g: If the interface is at x=0; two models: y = x ; y = x2;
both have the same value Y= 0 @ x=0;
• For y=f(x) and y=g(x), C1 exists at the interface between
two models (x=x0), if the value of the y (prediction) is the
same ; df(x0)/dx = dg(x0)/dx
• E.g: If the interface is at x=0; two models: y = x ; y = x2
do not have C1 continuity@ x=0 as: for y=x; dy/dx = 1 and
for y= x2 ; dy(x=0)/dx = 0.
Enforcing C0 Continuity: Example
• Continuity can be enforced by reducing parameters:
• Say f(x) - y = a0+a1 x and y = b0+b1 x must be C0
continuous at x = 1;
• Enforcing will lead to reduction in 1 constant
– a0+a1 = b0+b1; i.e. b0 = (a0+a1)-b1
– Choosing the following two function will enforce C0 continuity
– f(x) = a0+a1 x ; g(x)= (a0+a1)-b1 (x-1)
