close

Enter

Log in using OpenID

embedDownload
Applications of MATLAB/Simulink for
Process Dynamics and Control
Yinjie Tang, Assistant professor
EECE, Washington University
(This lecture was modified from slides provided by Professor Kirk
Dolan and Wei Liao at MSU and Venkat Subramanian at WashU)
3/17/2014
2
1. MATLAB Introduction
 Matlab is short for “MATrix LABoratory”
 High-performance technical computing environment
 Comprehensive math
 Graphic functions
 Powerful high-level language
 Simulink is a platform for multidomain simulation and model-based
design of dynamic systems
 Process control dynamics (ChE462 teaches design of feedback and
feedforward control using Simulink)
 Bioprocess Engineering (Solve multiple kinetic models)
Background
Structure of Matlab/Simulink
Matlab Platform
MATLAB based software…
COMSOL
Simulink
Toolboxes
M. File
Formerly FEMLAB
Curvefitting toolbox
A finite element analysis,
solver and simulation
Optimization toolbox
Differential equation toolbox
Other toolboxes……
Systems biology toolbox (new)
MATLAB fundamentals
• MATLAB can automatically handle rectangular arrays of data - one-dimensional
arrays are called vectors and multi-dimensional arrays are called matrices.
• Arrays are set off using square brackets [ ] in MATLAB
• Entries within a row are separated by spaces or commas
• Rows are separated by semicolons
•The transpose operator (apostrophe) can be used to flip an array over its own
diagonal. For example, if b is a row vector, b’ is a column vector.
Examples:
Q: What is a’=?
5
MATLAB fundamentals
You can use help check more commands:
logspace linspace
Clearing Commands
• When running a program many times, the command window may become
cluttered. Clear the command window with ―clc. (clear command).
• Good programming practice: At the beginning of the program, clear all variables:
―clear: removes all variables from workspace (“clear all” clears all objects in
workspace, plus resets all assumptions.
MATLAB fundamentals
Hold on and hold off
– hold on keep the current
data plotted and add the
results of any further plot
commands to the graph.
This continues until the
hold off command, which
clear the graph and start
over if another plotting
command is given. hold on
should be used after the
first plot in a series is made.
subplot(m, n, p)
– subplot splits the figure window into an
m×n array of small axes and makes the
pth one active. Note - the first subplot is
at the top left, then the numbering
continues across the row. This is different
from how elements are numbered within
a matrix!
3/17/2014
7
MATLAB fundamentals
There are two main kinds of M-file
– Script files: A script file is merely a set of MATLAB commands that are
saved on a file - when MATLAB runs a script file, it is as if you typed the
characters stored in the file on the command window. Scripts can be executed
either by typing their name (without the .m) in the command window, by
selecting the Debug, or Run (or Save and Run) command in the editing
window.
– Function files: Function files serve an entirely different purpose from script
files. Function files can accept input arguments from and return outputs to the
command window, but variables created and manipulated within the function
do not impact the command window. A function file can contain a single
function, but it can also contain a primary function and one or more
subfunctions. The primary function is whatever function is listed first in the
M-file - its function name should be the same as the file name. Subfunctions
are only accessible by the main function.
3/17/2014
8
MATLAB fundamentals
The general syntax for a function is:
function outvar = funcname(arglist)
statements %help comments
outvar = value;
where
– outvar: output variable name
– funcname: function name
– arglist: input argument list
– helpcomments: text to show with
help funcname
– statements: MATLAB commands
for the function
Script and function files Debug
MATLAB fundamentals (Input and Output)
Data can be read
from Excel files
3/17/2014
10
MATLAB fundamentals (Input and Output)
3/17/2014
11
MATLAB fundamentals (Input and Output)
More options: See also xlswrite, csvread, csvwrite, dlmread, dlmwrite, textscan.
3/17/2014
12
MATLAB fundamentals (programming)
Rational Operators
3/17/2014
Example
if I == J
A(I,J) = 2;
elseif abs(I-J) == 1
A(I,J) = -1;
else
A(I,J) = 0;
end
13
MATLAB fundamentals (programming, Loops)
3/17/2014
2. Matlab ODE45 function
(solve kinetic models)
• Ordinary Differential Equation (ODE): The dependent
variable is a function of only one independent variable.
• Partial Differential Equation (PDE): The dependent
variable is a function of more than one dependent variable.
• The order of a differential equation is the highest derivative
in the expression.
• A differential equation is linear if the unknown function and
its derivatives appear to the power one without products,
and nonlinear otherwise.
Initial-Value Problem
ODE equations
dy/dt=f(t, y)
y(0)=initial value; Find y(t)
Use Excel to Solve
Use MATLAB ode45
Stiff system
1. The solutions to some ODE problems
exhibit multiple time scales - for some
parts of the solution the variable changes
slowly, while for others there are abrupt
changes. Constant step-size algorithms
would have to apply a small stepsize to
the entire computation, wasting many
more calculations on regions of gradual
change.
2. MATLAB’s ode23 function uses
second- and third order functions to
solve the ODE. ode45 function uses
fourth- and fifth order RK functions.
Ode23 is better than ode45 for stiff
problems.
MATLAB has a number of built-in functions for solving stiff systems of ODEs,
including ode15s, ode23s, ode23t, and ode23tb
Solve high order ODEs
(Express an nth order ODE as a system of n first-order ODEs)
Convert the second-order ODE into a pair of first order ODEs
by defining
2. Matlab Application (b kinetic models)
k2
k1
ES
E+S
b) kinetic models
E+P
k-1
d [S ]
= − k1[ S ][ E ] + k−1[ ES ]
dt
d [E ]
= −k1 [S ][E ] + (k −1 + k 2 )[ES ]
dt
d [ ES ]
= k1[ S ][ E ] − ( k−1 + k2 )[ ES ]
dt
d [ P]
= k2 [ ES ]
dt
Matlab solution to
equations using ode45
Kinetic model 1: Enzyme kinetics)
%=======================================================
% dynamic enzyme kinetics: filename "Enzymekinetics.m"
%=======================================================
clear
clf
global k1 k2 k3; % define rate constant, k3 is k1-minus.
k1=100; k2=0.1; k3=100;
% initial conditions
S0=1; E0=1; ES0=0; P0=0; tend=100;
y0=[S0 E0 ES0 P0];
[t,y]=ode45('MMmodel',[0 tend], y0);
% results
S=y(:,1); E=y(:,2); ES=y(:,3); P=y(:,4);
%********************************************
% filename MMmodel.m, This function is used to
% calculate the MM model using ODE45
% Initial conditions: y0=[S0 E0 ES0 P0]
%*******************************************
function ydot=MMmodel(t,y)
global k1 k2 k3;
ydot(1)= -k1*y(1)*y(2)+k3*y(3); %substrate
ydot(2)= -k1*y(1)*y(2)+(k3+k2)*y(3); %enzyme
ydot(3)= k1*y(1)*y(2)-(k3+k2)*y(3); %enzyme
complex
ydot(4)= k2*y(3); % product
subplot(211)
ydot=ydot‘;
plot (t, S, 'k', t, P, 'y');
title('Substrate vs Products'), xlabel('time'), ylabel('amount'),
legend('substrate','product', 1);
subplot(212)
plot (t, E, 'k', t, ES, 'y');
title('Enzyme vs Enzyme complex'), xlabel('time'), ylabel('amount'),
legend('Enzyme','Enzyme complex', 1);
Enzyme kinetic model
3. Solve Partial Differential Equations
using MATLAB solver ‘pdepe’
Lian He
3/17/2014
21
Syntax and description
MATLAB command:
sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)


∂u  ∂u
∂x  ∂t
=
c  x, t , u , 
x
pdefun
−m
∂ m 
∂u   
∂u 
+
x
f
x
t
u
s
x
t
u
,
,
,
,
,
,

 

∂x  
∂x   
∂x 
icfun
u ( x, t 0 ) = u0 ( x )
bcfun
∂u 

0
p ( x, t , u ) + q ( x, t ) f  x, t , u ,  =
∂x 

Example 1
∂u ∂  ∂u 
π
=  
∂t ∂x  ∂x 
2
Initial condition
0≤x ≤1, t≥0
u ( x,0 ) = sin π x
Boundary conditions
u ( 0, t ) = 0
π e−t +
∂u
0
(1, t ) =
∂x
PDE function
∂u ∂  ∂u 
π
=  
∂t ∂x  ∂x 
2
∂u  ∂u − m ∂  m 
∂u   
∂u 

=
+
c  x, t , u , 
x
x
f
x
t
u
s
x
t
u
,
,
,
,
,
,




∂x  ∂t
∂x 
∂x   
∂x 


MATLAB commands:
function [c,f,s] = pdefun(x,t,u,DuDx)
c = pi^2;
f = DuDx;
s = 0;

Initial condition
u ( x, t 0 ) = u0 ( x )
u ( x,0 ) = sin π x
MATLAB commands:
function u0 = pdeic(x)
u0 = sin(pi*x);

Boundary conditions
u ( 0, t ) = 0
∂u 
u   0  

∂u 
p
x
,
t
,
u
+
q
x
,
t
f
x
,
t
,
u
,
0
( ) ( ) 
0
=
∂u
 −t  + 1  f  x, t , u, ∂x  =
−t
∂x 

0
π e + (1, t ) =

π e    
∂x
MATLAB commands:
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;

Call the solver ‘pdepe’

MATLAB commands:
m = 0;
xrange = linspace(0,1,20);
trange = linspace(0,2,20);
sol =
pdepe(m,@pdefun,@pdeic,@pdebc,xrange,trange);
xrange = linspace(0,1,10);
trange = linspace(0,2,10);
xrange = linspace(0,1,20);
trange = linspace(0,2,20);
Example 2
∂u1 ∂ 2u1
=
+ u1 (1 − u1 − u 2 )
2
∂t
∂x
∂u2 ∂ 2u2
=
+ u2 (1 − u1 − u 2 )
2
∂t
∂x
Initial
conditions:
Boudnary
conditions:
u1 ( 0, x ) = x 2
u2 (0, =
x) x(x − 2)
∂u1
=
u2 (t ,0) 0
( t ,0 ) 0,=
∂x
∂u2
=
u1 ( t ,1) 1,=
( t ,1) 0
∂x
PDE function
∂u1
=
∂t
∂u2
=
∂t
∂ 2u1
+ u1 (1 − u1 − u 2 )
2
∂x
∂ 2 u2
+ u2 (1 − u1 − u 2 )
2
∂x
1 ∂ u1  ∂  ∂ u1   u1 (1 − u1 − u 2 ) 
  + 


1 ⋅ ∂t =
  u2  ∂x  ∂x u2   u2 (1 − u1 − u 2 ) 
∂u  ∂u − m ∂  m 
∂u   
∂u 

=
c  x, t , u , 
x
x f  x, t , u ,   + s  x, t , u , 

∂x  ∂t
∂x 
∂x   
∂x 


Initial condition
u1 ( 0, x ) = x 2
u2 (0, =
x) x(x − 2)
 x2

u ( 0, x ) = 

 x(x − 2) 
Boundary conditions
∂u1
=
u2 (t ,0) 0
( t,0 ) 0,=
∂x
∂u2
=
u1 ( t ,1) 1,=
( t,1) 0
∂x
0  1 
0
u  +  0  f =
∂u 

 2  
p ( x, t , u ) + q ( x, t ) f  x, t , u ,  =
0
∂x 

u1 − 1 0 
0
0  + 1  f =

  
MATLAB program
function [c,f,s]=pdex2pde(x,t,u,DuDx)
c = [1;1];
f = [1;1].*DuDx;
s = [u(1)*(1-u(1)-u(2));
u(2)*(1-u(1)-u(2))];
sol=pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);
% ----------------------------------m=0;
x=linspace(0,1,20);
t=linspace(0,1,20);
u1=sol(:,:,1);
u2=sol(:,:,2);
subplot(2,1,1)
surf(x,t,u1)
xlabel('Distance x')
ylabel('Time t')
subplot(2,1,2)
surf(x,t,u2)
xlabel('Distance x')
ylabel('Time t')
function u0=pdex2ic(x)
u0 = [x^2;
x*(x-2)];
% -----------------------------------------function [pl,ql,pr,qr]=pdex2bc(xl,ul,xr,ur,t)
pl = [0;ul(2)];
ql = [1;0];
pr = [ur(1)-1;0];
qr = [0;1];
Results
PDE Toolbox can solve a general PDE
∂ 2u
∂ 2u
∂ 2u ∂u ∂u
A 2 + 2B
+C 2 + D + E + F =
0
∂x
∂x∂y ∂y
∂x ∂y
4. Parameter Estimation

Linear estimation:
The relationship between the dependent variables (y) and
independent variables (x) is linear
For example, y=a0+a1*x1+a2*x2

Nonlinear estimation:
The relationship between the dependent variables (y) and
independent variables (x) is nonlinear
For example, y=a0+a1*x+a2*x2+b0*xb1
Nonlinear Parameter Estimation: Monod Model
S mmolar/L
0
100
200
300
400
500
600
700
800
900
1000
1100
1200
1300
1400
1500
1600
1700
1800
1900
2000
V (mmolar/min)
0.0147
0.0488
0.1067
0.1420
0.2045
0.2094
0.2414
0.2161
0.2434
0.2188
0.2493
0.2546
0.2776
0.2563
0.2681
0.2718
0.2857
0.3132
0.2771
0.2738
0.2885
Vmax ⋅ S
V =
K +S
Vmax, and K are the unknown parameters
that need to be estimated
Curve Fitting Toolbox in MATLAB
Pros
1. User-friendly
2. Statistic analysis embedded
3. Multiple options available
Cons
1. Only one dependent variable allowed
2. Defined algebraic equations
3. Limited copies
Find more about Curve Fitting Toolbox in
http://www.mathworks.com/products/curvefitting/
Curve Fitting Toolbox in MATLAB
Input the data
Curve Fitting Toolbox in MATLAB
Choose Custom Equations
Fitting results
Define Custom Equations
Curve Fitting Toolbox in MATLAB
Choose different initial guess
Nonlinear Parameter Estimation: Monod Model
S mmolar/L
0
100
200
300
400
500
600
700
800
900
1000
1100
1200
1300
1400
1500
1600
1700
1800
1900
2000
V (mmolar/min)
0.0147
0.0488
0.1067
0.1420
0.2045
0.2094
0.2414
0.2161
0.2434
0.2188
0.2493
0.2546
0.2776
0.2563
0.2681
0.2718
0.2857
0.3132
0.2771
0.2738
0.2885
Vmax ⋅ S
V =
K +S
Vmax = 0.3463 [0.3176 0.3750]
K = 376.8 [265.2 488.4]
Example 2 (curve fitting toolbox)*
Example
 Creep compliance of a wheat protein film (determination of retardation
time and free dashpot viscosity in the Jefferys model)
J = J1 (1 − exp(
−t
λret
)) +
t
µ0
Note: creep is the tendency of a solid
material to slowly move or deform
permanently under the influence of stresses.
Where J is the strain, J1 is the retarded compliance (Pa-1); λret=µ1/G1 is retardation
time constant (s); µ0 is the free dashpot viscosity (Pa s); t is the time.
Example 2
Parameter estimation - Creep compliance Experimental Data
Note: Unlike brittle fracture,
creep deformation does not
occur suddenly upon the
application of stress. Instead,
strain accumulates as a result
of long-term stress. Creep is a
"time-dependent"
deformation.
Chewing gum
t, Time (s)
300
600
900
1200
1500
1800
2100
2400
2700
3000
3300
3600
3900
4200
4500
4800
5100
5400
5700
6000
6300
6600
6900
7200
J, Strain (Mpa^-1)
0.17
0.265
0.318
0.348
0.366
0.376
0.382
0.386
0.388
0.39
0.392
0.393
0.395
0.396
0.397
0.398
0.4
0.401
0.402
0.403
0.404
0.405
0.406
0.408
Creep compliance
1. Input data
3. Define Model
2. Select Data
Initial guesses and
boundaries are important
for a global solution
4. Initial Guess
Parameter estimation - Creep compliance
Results
Initial value (start point): a = 0.1, b = 1000, c = 100000
The results:
J1 = a = 0.38 Mpa-1
λret = b = 510.6 s
µ0 = c = 260800 Mpa s
J = 0.38(1 − exp(
−t
t
)) +
510.6
260800
The procedure:
http://www.swarthmore.edu/NatSci/echeeve1/Ref/MatlabCurveFit/MatlabCftool.html
Summary
Different ways to estimate parameters
 excel “solver”
A graphical user interface (curve fitting
toolbox)
 Matlab optimization tools (nlinfit or
fmincon)
3/17/2014
47
Introduction to Optimization
The maximizing or minimizing of a given function subject to
some type of constraints. Make your processes as effective as
possible.
In a typical chemical plant, there are many control variables for controlling the
process, such as maintaining a temperature, level, or flow.
Process Control and Optimization


Control has to do with
adjusting flow rates to
maintain the controlled
variables of the process at
specified set-points.
Optimization chooses the
values for key set-points
such that the process
operates at the “best”
economic conditions.
What is optimal operation
temperature?
Economic Objective Function


VB > VC, VA, or VAF
V is the chemical values.
At low T, little formation of B; At high T, too much of B reacts to
form C; Therefore, the exits an optimum reactor temperature, T*
Φ = Q C A VA + Q C B VB + Q CC VC − Q C A0 VAF
Graphical Solution of Optimum Reactor Temperature, T*
k1(T)
k2(T)
A  B  C
Use of “fmincon” in MATLAB for process optimization
ceq, and c describes the nonlinear
equalities/inequalities among parameters
A and b describes the linear inequalities
among parameters
Aeq and beq describes the linear equalities
among parameters
X = fmincon(fun, X0, A, B, Aeq, Beq, lb, ub) defines a set of lower and upper bounds on the design
variables, X, so that a solution is found in the range LB <= X <= UB. Use empty matrices for LB and
UB if no bounds exist.
[x fval] = fmincon(fun,x0,A,B,Aeq,Beq,lb,ub,nonlcon,options)
Estimated
parameters
Errors of predicted Function that you Initial guess of
and observed
want to minimize, the parameters
results
i.e. the error
Boundaries of
the parameters
Choose the right solver
The function NONLCON accepts X and returns the vectors
C and Ceq, representing the nonlinear inequalities and
equalities respectively.
Example 1
Find values of x that minimize
starting at the point x = [10; 10; 10] and subject to
the constraints
X = fmincon(fun, X0, A, B, Aeq, Beq, lb, ub) minimizes FUN subject to the
linear equalities Aeq*X = Beq as well as A*X <= B. (Set A=[] and B=[] if no
inequalities exist.)
Use fmincon for optimization
%file name example1
clc;close all;clear all;
int_guess=[10;10;10];% INITIAL GUESS FOR U'S
A=[-1 -2 -2;1 2 2]; % equality constraints
B=[0 72];
[u]=fmincon(@eg_1,int_guess,A,B) % CALL FOR OPTIMIZER
%file name eg_1.m
function [err]=eg_1(u)
x1=u(1);
x2=u(2);
x3=u(3);
err=-x1*x2*x3;
Example 2 (dynamic optimization)

Consider the batch reactor with following reaction
A
B
C
Find the temperature, at which the product B is maximum
Mathematical Representation of the system is as :
clc;close all;clear all;
warning off
int_guess=300;% INITIAL GUESS FOR U'S
LB=298; % LOWER BOUND OF U
UB=398;% UPPER BOUND ON U
[u,FVAL]=fmincon(@reactor_problem,int_guess,[],[],[],[],LB,UB,[]); %
CALL FOR OPTIMIZER PROBABLY NOT TO CHANGE BE USER
1
[err]=reactor_problem(u);
load data_for_system
plot(t,Y(:,1),'k')
hold on
plot(t,Y(:,2),'m')
legend('opt c_A','opt c_B')
u
opt c A
opt c B
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
T=335.3244
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Call ODE to solve the model
ode45 (@fun, slot, init, [opt], [par1, par2, ...])
Maximize B
function [err]=reactor_problem(u)
A=[1 0]; %initial condition
[t,Y]=ode15s(@reactor_ODE,[0 1],A,[],u); %tspan=1
err=-Y(end,2); %make the maximum B
save data_for_system t Y %just for plotting
ODE calculation
function [YPRIME]=reactor_ODE(t,x,u)
YPRIME=zeros(2,1);
T=u;
k1=4000*exp(-2500/T);
k2=620000*exp(-5000/T);
YPRIME(1)=-k1*(x(1))^2;
YPRIME(2)=(k1*x(1)^2)-k2*x(2);
Part 2: Inverse problem (from experimental
data to model construction)
Parameter fitting for dynamic model
An example for parameter fitting of
dynamic bioprocess model
A) Nonlinear Parameter Estimation via “fmincon”
in MATLAB
References (e.g.
literature)
Mathematic
Model
References (e.g.
literature)
Mathematic
Model
v.s.
Known
Parameters
Unknown
Parameters
Simulated
Results
“Parameter
estimation”
ε
v.s.
“Forward problem”
Experimental
Results
Simulated
Results
“Inverse problem”
Experimental
Results
Nonlinear Parameter Estimation
1. Parameter estimation for linear systems
y = A·x + b
y, x are the vector of dependent and independent variables,
A, b are parameters to be estimated
“regress” command in MATLAB (for linear model)
2. Parameter estimation for nonlinear systems
y = f (parameters, x)
f is the nonlinear function of the estimated parameters
For example, y=β1+β2·sin(β3· t)
cftool box for more complicated model equations
3. Parameter estimation for dynamic ode models (“fmincon” , “nlinfit”, etc.
in MATLAB)
Example 1: parameter estimation in bioengineering
• Bioreactor model
dX
= µ⋅X
dt
dP
= YP / X ⋅ µ ⋅ X
dt
1
dS
=−
⋅µ⋅ X − µ⋅ X
dt
YX / S
µ=
µmax ⋅ S
KS + S
Five parameters to be estimated:
YP/X, YX/S, µmax, Ks, and S(0)
Initial conditions:
X(0) = 0.05 g/L; P(0) = 0 g/L;
Experimental data
t
0.0
3.0
6.0
9.0
12.0
15.0
18.0
21.0
24.0
27.1
X
0.51
0.15
0.42
0.56
0.33
0.09
1.01
1.60
2.62
3.73
P
0.28
0.47
0.21
0.61
0.12
0.09
0.16
0.34
0.52
1.13
S
14.86
14.93
14.30
14.00
13.80
12.00
11.11
6.00
3.50
0.70
We will use “fmincon” to find the
parameters
The Error Function
function [error,ypred] = reactor_model(beta)
global yobs
global t
Input the observed data
%reactor model
S0=beta(5);
y0=[0.05 0 S0];
tspan=t; %we want y at every t
[t,y]=ode45(@ff,tspan,y0);
Settings to run ODEs
function dy = ff(t,y) %function that computes the dydt
umax=beta(1);
Ks=beta(2);
Ypx=beta(3);
Yxs=beta(4);
X=y(1);
P=y(2);
S=y(3);
u=umax.*S./(Ks+S);
dy(1)=u.*X;
dy(2)=Ypx.*u.*X;
dy(3)=-1/Yxs.*u.*X-u.*X;
dy=dy';
end
y1=y(:,1); y2=y(:,2); y3=y(:,3);
ypred=[y1; y2; y3];
error=sum((ypred-yobs).^2);
end
The ODEs that describe
bioreactor behaviors
Calculate the error
The fmincon Function
%Use fmincon to find the parameters in reactor model
clear
clc
data =xlsread('reactor_model_data.xlsx');
t=data(:,1);
X=data(:,2);
P=data(:,3);
S=data(:,4);
x=t;
Read data from Excel
global yobs
global t
yobs=[X;P;S];
umax=0.1;
Ks=10;
Ypx=0.11;
Yxs=0.45;
S0=14.86;
beta0(1)=umax; %initial guess
beta0(2)=Ks; %initial guess
beta0(3)=Ypx; %initial guess
beta0(4)=Yxs; %initial guess
beta0(5)=S0; %initial guess
Set up the initial guess of
parameters
The fmincon Function (continued)
%fmincon
fun=@reactor_model;
A=[];
b=[];
Aeq=[];
beq=[];
lb=zeros(5,1);
ub=lb+20;
nonlcon=[];
x0=beta0;%lb+rand.*(ub-lb);
options=optimset('Algorithm','interior-point');
Settings to run fmincon
Q: How to find the global solution?
[param, fval] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
%get the predicted resutls
[error,ypred]=reactor_model(param);
n=length(t);
Xpred=ypred(1:n);
Ppred=ypred(n+1:2*n);
Spred=ypred(2*n+1:3*n);
figure(1)
set(gca, 'fontsize',14,'fontweight','bold');
hold on
plot(x,Xpred,'-r','linewidth',2.5);
plot(x,Ppred,'-g','linewidth',2.5)
plot(x,Spred,'-b','linewidth',2.5)
plot(x,X,'r+','markersize',10);
plot(x,P,'go','markersize',10)
plot(x,S,'bx','markersize',10)
xlabel('time')
ylabel('y')
legend('X g/L','P g/L','S g/L')
A: Randomly perturb the initial guess
Run the model again to
get the predicted results
Plot the predicted v.s. observed
results
Parameters Estimated via “fmincon”
param =
0.2248 3.8128 0.2606 0.3029 14.9987
fval =
4.3559
µmax = 0.2248 (1/h)
Ks = 3.8128 (g/L)
Ypx = 0.2606 (g/g)
Yxs = 0.3029 (g/g)
S0 = 14.9987 (g/L)
Estimating confidence intervals of parameters via
“fmincon”
Which one makes sense? Numbers in the brackets are the 95% confidence interval of µmax
µmax = 0.2248 [0.1700 0.2797] (1/h)
Or
µmax = 0.2248 [-0.1064 0.2874] (1/h)
Estimation of µmax is reliable
µmax cannot be accurately estimated
How to estimate the confidence intervals of parameters?
• “bootstrap” method: randomly resample the experimental observations; and
for each new set of sampled experimental observations, use fmincon to find a
new set of parameters.
• “Monte Carlo” method: randomly perturb the experimental observations within
the measurement errors (i.e. standard derivation); and for each new set of
perturbed experimental observations, use fmincon to find a new set of parameters.
B: Parameter fitting using nlinfit
Unlike forward problems, inverse problems require experimental data, and an
iterative solution. Because inverse problems require solving the forward
problem numerous time, the ode45 solver will be nested within a nonlinear
regression routine called “nlinfit.”
The syntax is: [param, r, J, COVB, mse] = nlinfit(X, y, fun, beta0);
returns the fitted coefficients param, the residuals r, the Jacobian J of function
fun, the estimated covariance matrix COVB for the fitted coefficients, and an
estimate MSE of the variance of the error term.
X is a matrix of n rows of the independent variable
y is n-by-1 vector of the observed data
fun is a function handle to a separate m-file to a function of this form:
yhat = fun(b,X)
where yhat is an n-by-1 vector of the predicted responses, and b is a vector of
the parameter values. beta0 is the initial guesses of the parameters.
Nonlinear Parameter Estimation: Bioprocess
µ ⋅S
dX
= µ ⋅ X = max
dt
KS + S
t (h)
0.0
3.0
6.0
9.0
12.0
15.0
18.0
21.0
24.0
27.1
X (g/L) P (g/L)
0.51
0.28
0.15
0.47
0.42
0.21
0.56
0.61
0.33
0.12
0.09
0.09
1.01
0.16
1.60
0.34
2.62
0.52
3.73
1.13
S (g/L)
14.86
14.93
14.30
14.00
13.80
12.00
11.11
6.00
3.50
0.70
dP
= YP / X ⋅ µ ⋅ X
dt
1
dS
=−
⋅µ ⋅ X − µ ⋅ X
dt
YX / S
Q: What is the problem for the model?
X: biomass,
P: product,
S: substrate
Initial conditions:
X(0)= 0.05 g/L
P(0)= 0 g/L
5 unknown parameters:
1. μmax (hr-1)
2. KS (g/L)
3. YP/X (g/g)
4. YX/S (g/g)
5. S(0) (g/L)
“nlinfit” in MATLAB
Pros
1. Multiple variables allowed
2. Customer defined equations
3. No limitations in copies
Cons
1. Programming by yourself
2. No direct statistic analysis
Syntax
beta = nlinfit(X,y,fun,beta0)
[beta,r] = nlinfit(X,y,fun,beta0)
[beta,r,J] = nlinfit(X,y,fun,beta0)
[beta,r,J,COVB] = nlinfit(X,y,fun,beta0)
[beta,r,J,COVB,mse] = nlinfit(X,y,fun,beta0)
[beta,...] = nlinfit(X,y,fun,beta0,options)
Output
beta: estimated
parameters
r: residual of the fit
(i.e. difference between
simulated data and real
data)
Find more about nlinfit in
http://www.mathworks.com/help/toolbox/stats/nlinfit.html
or searching “nlinfit” in product help of MATLAB
Input
X: independent variables
y: dependent variables
fun: Bioprocess equations
beta0: initial guess of
unknown Parameters
1/--pages
Report inappropriate content