Transcript Part B

Inverse problem (from experimental data to
model construction)
Part
2
Parameter fitting for ODEs using fmincon function
Example by Xueyang Feng: Nov 16th
X = FMINCON(FUN,X0,A,B,Aeq,Beq) minimizes FUN subject to the linear
equalities Aeq*X = Beq as well as A*X <= B. (Set A=[] and B=[] if no
inequalities exist.)
1
4/8/2015
d) Parameter fitting using nlinfit
Inverse problem
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.
120
100
80
60
Based on experimental
data, find y0 and k.
y
Example 1: Model fitting ODE
equation
dy
ky
dt
ypred
yobs
40
20
data =xlsread('exp_data.xls'); %read data from excel 0
yo=100; k=0.6;
-20
0
1
2
3
4
5
6
7
beta0(1)=yo;
time (min)
beta0(2)=k;
x=data(:,1); yobs=data(:,2);
[param,resids,J,COVB,mse] = nlinfit(x,yobs,'forderinv',beta0);
rmse=sqrt(mse); %root mean square error = SS/(n-p)
%R is the correlation matrix for the parameters, sigma is the standard error vector
[R,sigma]=corrcov(COVB);
%confidence intervals for parameters
ci=nlparci(param,resids,J);
%computed Cpredicted by solving ode45 once with the estimated parameters
ypred=forderinv(param,x);
%mean of the residuals
meanr=mean(resids);
8
9
10
figure
hold on
h1(1)=plot(x,ypred,'-','linewidth',3); %predicted y values
h1(2)=plot(x,yobs,'square', 'Markerfacecolor', 'r');
legend(h1,'ypred','yobs')
xlabel('time (min)')
ylabel('y')
%residual scatter plot
figure
hold on
plot(x, resids, 'square','Markerfacecolor', 'b');
YLine = [0 0];
XLine = [0 max(x)];
plot (XLine, YLine,'R'); %plot a straight red line at zero
ylabel('Observed y - Predicted y')
xlabel('time (min)‘)
Function with ode45
function y = forderinv(param,t)
%first-order reaction equation
tspan=t; %we want y at every t
[t,y]=ode45(@ff, tspan, param(1));
%param(1) is y(0)
function dy = ff(t, y) %function that
computes the dydt
dy(1)= -param(2)*y(1);
end
end
Raw data
120
ypred
yobs
100
80
One ODE
y
60
40
20
0
-20
0
1
2
3
4
5
6
time (min)
7
8
9
10
15
Observed y - Predicted y
10
5
0
-5
-10
-15
0
1
2
3
4
5
6
time (min)
7
8
9
10
Time
0
0.10101
0.20202
0.909091
1.010101
1.111111
1.212121
1.313131
1.616162
1.717172
1.818182
1.919192
2.020202
2.121212
2.525253
2.626263
2.727273
2.828283
2.929293
3.232323
3.333333
3.434343
3.535354
3.939394
4.040404
4.343434
4.444444
4.545455
4.646465
4.747475
4.848485
4.949495
y
107.2637
102.9394
92.71275
75.15614
70.4741
69.83605
57.00548
60.03961
54.01795
56.79646
53.81866
49.67512
42.0459
40.42633
41.06735
39.48581
34.28254
30.3689
31.69877
25.20435
27.20177
19.71342
26.31708
23.20818
15.5449
19.08923
16.08262
19.26572
23.76179
11.85451
7.628606
7.998508
File name: exp_data.xls
Time
5.050505
5.454545
5.555556
5.656566
5.757576
5.858586
5.959596
6.060606
6.161616
6.262626
6.666667
6.767677
6.868687
6.969697
7.070707
417395
7.575758
7.676768
7.777778
8.181818
8.282828
8.383838
8.484848
8.888889
8.989899
9.090909
9.191919
9.292929
9.393939
9.494949
9.59596
9.69697
10
y
8.541471
12.46344
6.944696
15.90549
5.717727
9.637526
4.531673
5.446719
7.203193
7.02317
16.22408
5.286706
11.7421
-4.34103
9.103846
9.998128
6.737791
7.460489
4.723043
8.394924
-0.45654
0.910225
7.662295
-0.04558
2.102022
1.454657
4.797709
9.16223
-5.94742
12.27148
5.956534
3.329466
Three ODES for batch fermentation
S
15
10
y
/s
X
P
5
Yp/s=Yx/s*Yp/x
0
0
5
10
15
time (min)
20
25
Three ODEs parameter fitting
File name: HW217.xls
t
X
P
S
0
0.05
0
15
2.7
0.082555
0.019533
14.86978
5.4
0.136259
0.051755
14.65496
8.1
0.224763
0.104858
14.30095
14
0.667768
0.370661
12.52893
20.5
2.151108
1.260665
6.595568
24.1
3.652464
2.161478
0.590144
Fitting four parameters
Script
clear all %clear all variables
global y0
data =xlsread('HW217.xls');
%initial conditions
y0=[0.05 0 15];
%initial guess
Umax=0.1; Ks=10; Ypx=0.11; Yxs=0.45;
beta0(1)=Umax; beta0(2)=Ks; beta0(3)=Ypx; beta0(4)=Yxs;
%Measured data
x=data(:,1);
yobsX=data(:,2);
yobsP=data(:,3);
yobsS=data(:,4);
yobs=[yobsX; yobsP; yobsS];
%nlinfit returns parameters, residuals, Jacobian (sensitivity coefficient matrix),
%covariance matrix, and mean square error. ode45 is solved many times
%iteratively
[param,resids,J,COVB,mse] = nlinfit(x, yobs,'forderinv2', beta0);
Script
rmse=sqrt(mse); %root mean square error = SS/(n-p)
n=size(x); nn=n(1);
%confidence intervals for parameters
ci=nlparci(param,resids,J);
%computed Cpredicted by solving ode45 once with the estimated
parameters
ypred=forderinv2(param,x);
ypredX=ypred(1:nn);
ypredP=ypred(nn+1:2*nn);
ypredS=ypred(2*nn+1:3*nn);
%mean of the residuals
meanr=mean(resids);
Script
figure
hold on
plot(x, ypredX, x, ypredP, x, ypredS); %predicted y values
plot(x, yobsX, 'r+', x, yobsP, 'ro', x, yobsS, 'rx');
xlabel('time (min)')
ylabel('y')
%residual scatter plot
x3=[x; x; x];
figure
hold on
plot(x3, resids, 'square', 'Markerfacecolor', 'b');
YLine = [0 0];
XLine = [0 max(x)];
plot (XLine, YLine,'R'); %plot a straight red line at zero
ylabel('Observed y - Predicted y')
xlabel('time (min)‘)
Function and sub-function
function y = forderinv2(param,t)
%first-order reaction equation
global y0;
tspan=t; %we want y at every t
[t,y]=ode45(@ff,tspan,y0); %param(1) is y(0)
function dy = ff(t,y) %function that computes the dydt
dy(1)= param(1)*y(3)/(param(2)+y(3))*y(1); %biomass
dy(2)= param(1)*param(3)*y(3)/(param(2)+y(3))*y(1); %product
dy(3)= -1/param(4)*dy(1)-1/(param(3)*param(4))*dy(2); %substrate
dy=dy';
end
% after the ode45, rearrange the n-by-3 y matrix into a 3n-by-1 matrix
% and send that back to nlinfit.
y1=y(:,1); y2=y(:,2); y3=y(:,3);
y=[y1; y2; y3];
end
3. Simulink
Basic elements
 Blocks
 Blocks: generate, modify, combine, and display signals
 Typical blocks
 Continuous: Linear, continuous-time system elements (integrators, transfer
functions, state-space models, etc.)
 Discrete: Linear, discrete-time system elements (integrators, transfer functions,
state-space models, etc.)
 Functions & Tables: User-defined functions and tables for interpolating
function values
 Math: Mathematical operators (sum, gain, dot product, etc.)
 Nonlinear: Nonlinear operators (coulomb/viscous friction, switches, relays,
etc.)
 Signals: Blocks for controlling/monitoring signals
 Sinks: Used to output or display signals (displays, scopes, graphs, etc.)
 Sources: Used to generate various signals (step, ramp, sinusoidal, etc.)
 Lines
 transfer signals from one block to another
3. Simulink
Simulink
Type “simulink” in
the command window
The graphic interfaces
Matlab
Tutorial example
x = sin (t)
simout
4/8/2015
Simulink example 2
Creep compliance of a wheat protein film
Using formaldehyde cross-linker
 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
Where J is the strain, J1 is the retarded compliance (Pa-1); λret=µ1/G1 is
retardation time (s); µ0 is the free dashpot viscosity (Pa s); t is the time.
 The recovery of the compliance is following the equation (t >t1):
J  J1 exp(
t  t1
ret
)
Where t1 is the time the stress was released.
4. Simulink examples
Creep compliance of a wheat protein film
Simulink model
Parameters: J1 = a = 0.38 Mpa-1; λret = b = 510.6 s; µ0 = c = 260800 Mpa s
J  J1 (1  exp(
t
ret
)) 
J  J1 exp(
t
0
t  t1
ret
)
4. Simulink examples
Creep compliance of a wheat protein film
Simulation result
4. Simulink example
Fermentation system
Agitation
Regular fermenter
Aeration
Ethanol fermenter
Heater
Heater
Products
Ethanol
Biomass
Biomass
4. Simulink examples
Ethanol production kinetics
  f ( s, p )
Growth kinetics

maxs
KS  s
(1 
p
pmax
)
 s
dx
p 
 rx  x  max (1 
) x
dt
KS  s
pmax
Substance consumption
dCS
max s

p 
 rS  [
 mS ]x  [
(1 
)  mS ]x
dt
YXS
YXS ( K S  s)
pmax
Product production
 s
dCP
p 
 rP  [YPX   mP ]x  [YPX max (1 
)  mP ]x
dt
KS  s
pmax
Where: η is the toxic power, Pmax is the maximum product concentration at which the growth is completely inhibited
4. Simulink examples
a. Using blocks to construct the model
Using S-function to construct the model (you may also writing the
programs for simulink function)
Special help session
November 28 (Monday, 6~9pm)
I will be in Sever 201 computer Lab
to answer your modeling questions.
4/8/2015
21
Introduction to Optimization
The maximizing or minimizing of a given function subject to some
type of constraints.
Make your processes as effective as possible
c) Introduction to Optimization
Example 1 (fminsearch)
To find the minimum of a multidimensional function in [-1.2, 1]:
f ( x)  100( x2  x12 ) 2  (1  x1 ) 2
The plot of the function can be generated by the following code:
[x, y]=meshgrid ([-2:0.1:2]);
z=100*(y-x.^2).^2+(1-x).^2;
mesh(x, y, z);
grid on
The function of ex2c is written as follow:
function f=ex2c(x)
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
The function of ex2cMinimum is written as follow:
function ex2cMinimum
fminsearch(@ex2c, [-1.2,1]);
clear
clc
[x, y]=meshgrid ([-2:0.1:2]);
z=100*(y-x.^2).^2+(1-x).^2;
mesh(x, y, z);
grid on
a=fminsearch(@ex2c, [-1.2,1]);
x=a(1)
y=a(2)
z=100*(y-x.^2).^2+(1-x).^2
%===========================================
% Michaelis & Menten Model Fit.
% File name "SimpleMMplot.m".
%===========================================
Example 2
Using fmin for curve fitting
clear
clf
S 
K m  S 
V  Vmax
global S V;
% experiemntal data
S=[0 1 2 3 4 5 6 7 8];
V=[0 0.08 0.15 0.18 0.2 0.21 0.215 0.216 0.216];
v0=[1; 1];
a=fminsearch('two_var',v0); %least square errors
Km=a(1);
Vmax=a(2);
Vmodel=Vmax*S./(Km+S);
plot(S,V,'ro', S, Vmodel)
function sumsqe= two_var(aa)
xlabel('Concentrations')
global S V;
ylabel('rate')
Km =aa(1); Vmax =aa(2);
legend('Data','Predict')
Vmodel=Vmax*S./(Km+S);
err=V-Vmodel; err2=err.*err;
sumsqe=sum(err2);
return
Example 3

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
FMINCON attempts to solve problems of the form:
min F(X) subject to: A*X <= B, Aeq*X = Beq (linear constraints)
C(X) <= 0, Ceq(X) = 0 (nonlinear constraints)
LB <= X <= UB
(bounds)
%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 4 (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
T=335.3244
0.1
0
FVAL = -0.6058
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
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);
Using MATLAB for solving more
complicated dynamic models
(Optional)
4/8/2015
30
For example, solve complicated
dynamic model in the bioprocess
Many biological systems are dynamic and
heterogeneous! Therefore, structured and segregated
models have to be used.
Definitions
• 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.
•Boundary-value problems are those where space
conditions are not known; Initial-value problems are those
dynamic condition are not known.
Boundary-value problems
• Dirichlet boundary conditions are those where a fixed value of a variable is
known at a particular location.
• Neumann boundary conditions are those where a derivative is known at a
particular location.
Finite-Difference Methods are often used to solve boundary-value problems.
In these techniques, finite differences are substituted for the derivatives in the
original equation, transforming a linear differential equation into a set of
simultaneous algebraic equations.
Final
equation
Marine pollution by oil compounds
Oil transport and degradation in
the porous sediment
4/8/2015
35
How to solve PDEs using the finite difference method?
Dm
rmaxC b
C b
 2C b 1 C b  2C b

[


]
2
2
t
(1   )  s K P   r
r r
Z
K s [(1   )  s K P   ]  C b
Ci , j

Z
The change of
Contaminant
concentrations in
the porous
environment
a: one injection unit
 2 Ci , j
Z
2
 2 Ci , j
r 2
Ci , j
r
Find the time and space dependent change
of Ci,j
Ci , j 1  Ci , j 1
2hZ



Ci , j 1  2Ci , j  Ci , j 1
hZ2
Ci 1, j  2Ci , j  Ci 1, j
hr2
Ci 1, j  Ci 1, j
2hr
Ci , j
 f (r , Z , Ci , j )
t
function ydot=rhs2d2(t, y)
global dx D ii jj Km Vmax
for j=1:jj
for i=1:ii
vv=vv+1;
yy(i,j) = y(vv);
end
end
Sample codes for PDE problem.
for j=2:jj-1
for i=2:ii-1
R=-Vmax*yy(i,j)/(Km+yy(i,j));
Dffu1=(yy(i+1,j)-2*yy(i,j)+yy(i-1,j))/(dx*dx);
Dffu2=(1/(dx*(i-1)))*(yy(i+1,j)-yy(i-1,j))/2/dx;
Dffu3=(yy(i,j+1)-2*yy(i,j)+yy(i,j-1))/(dx*dx);
dot(i,j) = D*(Dffu1+Dffu2+Dffu3)+R;
bb = ii*(j-1)+i;
ydot(bb) = dot(i,j);
end
end
ydot=ydot';
4/8/2015
You can change
PDE problem to
many mini ODE
problem, then you
solve them at the
same time.
Solving Reaction Engineering Problems
with MATLAB
By Lian He
Solve PDE dynamic equations and
nonlinear equilibrium equations
Example 1

Porous Catalyst
The spherical catalysts are exposed to reactants at the
initial time point. Given the information about the
diffusion coefficient De, reaction rate, etc., we want to
know the reactant concentration profile. (We assume that
the reactant concentration in the catalyst is only a function
of time and distance to the center for simplicity. )
R
Solution
• PDE

Dimensionless form
Initial condition
u(0, x)=0
Thiele Module
Boundary conditions
u(t , 1)=Ssurf
MATLAB tool: pdepe

Syntax
sol = pdepe(m, @pdefun, @icfun, @bcfun, xmesh, tspan, options)
function [pl,[c,
ql, pr,
qr]=pdebound(xl,
function
f, s]
= pdefun(x,ul,
t, xr,
u, ur, t)
pr=ur-1;
DuDx)
qr=0;
function u0=icfun(x)
c=1;
pl=0;
u0=0;
ql=1;
f=φ^(-2)*DuDx;
end
end
s=-u/(1+u/beta);
end
Results
Main script
clc;
clear all;
x=0:0.02:1; %radius range
t=0:0.04:8; %timespan
u=pdepe(2,@pdefun,@icfun,@bcfun,x,t);%m=2
%Show the concentration profile as a function of time
for k=1:length(t)
plot(x,u(k,:),'linewidth',2);
xlabel('distance from the center');
ylabel('substrate concentration');
pause(.2)
end
figure(2)
%3D plot
mesh(x,t,u);
xlabel('x');
ylabel('t');
zlabel('u');
axis([0 1 0 8 0 1]);
Three Functions
function [pl, ql, pr, qr]=bcfun(xl, ul, xr, ur,
t)
pr=ur-1;
qr=0;
pl=0;
ql=1;
end
function u0=icfun(x)
u0=0;
end
function[c,f,s]=pdefun(x,t,u,DuDx)
De=3*10^-4;
Ss=1*10^-4;
Km=0.1*Ss;
R=0.001;
umax=0.5;
Mt=R*sqrt(umax/Km/De);
beta=Km/Ss;
c=1;
f=Mt^(-2)*DuDx;
s=-u/(1+u/beta);
end
Example 2

Solution
#
Reaction
Molar Extent
1
C+2H2O=CO2+2H2
X1
2
C+CO2=2CO
X2
3
C+2H2=CH4
X3
Species
nj0
nj
H2O
1
1-2X1
CO2
0
X1-X2
CO
0
2X2
CH4
0
X3
H2
0
2(X1-X3)
sum
1
1+X1+X2-X3
Nonlinear equations
You can use MATLAB tool-fsolve to solve the problem
clear all;
clc;
%initial guess of molar extent
x0=[0.3 0.3 0.3];
%set the intial values
T=600;
%creat different vectors for recording temperature, molar extent and molar
%fraction
t=zeros(101,1); X=zeros(101,3); yH2O=zeros(101,1); yCO2=zeros(101,1);
yCO=zeros(101,1); yCH4=zeros(101,1); yH2=zeros(101,1);
for i=1:101
%find the solution
[x,fval] = fsolve(@nlfunc,x0,[],T);
%record the results
t(i)=T; X(i,:)=x;
%caculate the molar fraction of each species
yH2O(i)=(1-2*x(1))/(1+x(1)+x(2)-x(3)); yCO2(i)=(x(1)-x(2))/(1+x(1)+x(2)-x(3));
yCO(i)=2*x(2)/(1+x(1)+x(2)-x(3)); yCH4(i)=x(3)/(1+x(1)+x(2)-x(3));
yH2(i)=2*(x(1)-x(3))/(1+x(1)+x(2)-x(3));
%change values for the next loops
T=T+10;
end
%plotting
subplot(2,2,[1 3]);
plot(t,X(:,1),'o',t,X(:,2),'*',t,X(:,3),'.');
legend('C+H_2O=>CO_2+2H_2','C+CO_2=>2CO','C+2H_2=>CH_4');
xlabel('Temperature (K)');
ylabel('molar extent'); axis([600 1600 0 0.6]);
subplot(2,2,[2 4]);
plot(t,yH2O,'o',t,yCO2,'-o',t,yCO,'.',t,yCH4,'*',t,yH2,'^');
legend('H_2O','CO_2','CO','CH_4','H_2');
xlabel('Temperature (K)');
ylabel('molar fraction');
axis([600 1600 0 0.6]);
Main script
Fsolve function to solve nonlinear equation
function F = nlfunc(x,T)
%coefficient values
K1=9.73*10^(-12)*exp(-10835/T+36.36);
K2=9.92*10^(-22)*exp(-20740/T+69.60);
K3=8.00*10^8*exp(8973/T-30.11);
%three column vectors
F=zeros(1,3);
%functions
F(1)=4*(x(1)-x(3))^2*(x(1)-x(2))-K1*(1-2*x(1))^2*(1+x(1)+x(2)-x(3));
F(2)=4*x(2)^2-K2*(x(1)-x(2))*(1+x(1)+x(2)-x(3));
F(3)=x(3)*(1+x(1)+x(2)-x(3))-4*K3*(x(1)-x(3))^2;
end
Results