Matlab Optimization

1
Matlab Optimization

• Optimization toolbox
• Solution of linear programs
• Metabolic flux balance analysis example
• Solution of nonlinear programs
• Batch fermentation example

2
Matlab Optimization Toolbox
Note fsolve is also part of the optimization
toolbox
3
Linear Programming (LP)

• Optimization of a linear objective function with
  linear equality and/or inequality constraints
• Standard LP form
• Matrix A must have more columns than rows
  (under-determined problem)
• Common solvers CPLEX, MOSEK, GLPK
• Further information

x vector of variables to be determined
(decision variables) A matrix of known
coefficients b vector of known coefficients c
vector of weights cTx scalar objective
function linear combination of the decision
variables
http//www-unix.mcs.anl.gov/otc/Guide/faq/linear-p
rogramming-faq.html
4
Matlab LP Solver linprog

• Solves linear programming (LP) problems of the
  form
• Syntax
• x linprog(f,A,b,Aeq,beq,lb,ub)
• Set Aand bif no inequality constraints
  exist
• Set Aeqand beqif no equality constraints
  exist
• Replace f with -f to find the maximum
• Defaults to a large-scale interior point method
  with options for a medium-scale simplex method
  variation or the simplex method
• See help linprog for additional details and
  options

where f, x, b, beq, lb, and ub are vectors and A
and Aeq are matrices
5
Metabolic Network Model

• Intracellular reaction pathways describing carbon
  metabolism
• Consumption of carbon energy sources (e.g.
  glucose)
• Conversion of carbon sources to biomass
  precursors (cell growth)
• Secretion of byproducts (e.g. ethanol)
• Each node corresponds to a metabolite
• Each path (line) corresponds to a reaction
• Stoichiometric matrix, A
• Row for each intracellular species (m rows)
• Column for each reaction (n columns)
• The entry at the ith row and jth column (ai,j)
  corresponds to the stoichiometric coefficient of
  species i participating in reaction j
• Av 0, stoichiometric balance on the metabolites
  where v is the vector of reaction fluxes
• More reactions (unknowns) than species
  (equations)
• Solution requires either enough measurements for
  the system to become square (n-m measurements) or
  optimization

6
Flux Balance Analysis (FBA)

• Linear programming (optimization) approach for
  resolving an under-determined metabolic network
  model
• Objective function based on an assumed cellular
  objective such as maximization of growth
• LP formulation
• Growth rate, m, represented as a linear
  combination of intracellular fluxes of the
  biomass precursors
• Flux bounds represent physiochemical or
  thermodynamic constraints on the reaction fluxes
• Extracellular conditions place limits on fluxes
  (e.g. oxygen availability)
• Thermodynamics constrain the direction a reaction
  may proceed reversible or irreversible
• The solution is the set of fluxes that maximizes
  cellular growth while satisfying the bounds and
  stoichiometric constraints

7
Flux Balance Analysis Example

• Yeast metabolic network model from HW 2
• Slightly modified to improve suitability for Flux
  Balance Analysis (FBA)
• 19x22 stoichiometric matrix
• Under-determined with 3 degrees of freedom
• Use FBA to determine solution corresponding to
  optimal cell growth

• Solve the LP
• gtgt v linprog(-w,,,A,b,vb(,1),vb(,2))
• Optimization terminated.
• View predictions for growth, oxygen uptake, and
  ethanol secretion
• gtgt mu w'v, vo2 v(io), ve v(ie)
• mu
• 101.9302
• vo2
• 108.3712
• ve
• 2.4147e-014
• All calculated values relative to a fixed glucose
  uptake rate normalized to 100

• Download the stoichiometric matrix to the Matlab
  working directory and load into Matlab
• gtgt load A.txt
• Specify the indices of key fluxes glucose,
  ethanol, oxygen, and biomass
• gtgt ig 22 ie 20
• gtgt io 19 imu 17
• Av 0
• gtgt m n size(A)
• gtgt b zeros(m,1)
• Objective function
• gtgt w zeros(n,1)
• gtgt w(imu) 1
• Specify flux bounds (all fluxes irreversible,
  glucose uptake fixed)
• gtgt vb zeros(n,1) Infones(n,1)
• gtgt vb(ig,) 100 100

8
FBA Example cont.

• Determine sensitivity of model predictions to the
  oxygen uptake rate to assess the tradeoff
  between achievable ethanol yields and cellular
  growth
• Create a vector of oxygen uptake rates to be
  considered
• gtgt vo 11125
• Implement a for loop to iterate over each entry
  in the oxygen uptake vector (vo). For each
  iteration (inside the loop), update the upper
  bound on oxygen uptake, solve the LP, and store
  the solution (mu, ve)
• gtgt for i1length(vo)
• vb(io,2) vo(i)
• v linprog(-w,,,A,b,vb(,1),vb(,2))
• mu(i) w'v
• ve(i) v(ie)
• end
• Plot the results
• gtgt plot(vo,mu,vo,ve)
• gtgt xlabel('Oxygen Flux')
• gtgt legend('Growth Rate','Ethanol Flux)
• Notice the tradeoff between cell growth and
  ethanol production. Highest ethanol
  productivity is achieved in batch fermentation by
  initially operating aerobically to rapidly
  increase cell density then switching to anaerobic
  conditions to produce ethanol.

A fixed oxygen uptake rate (lower bound equal to
upper bound) was not specified to avoid forcing
the cell to take up more oxygen than necessary
(cell should only use as much oxygen as it needs).
9
Nonlinear Programming (NLP)

• Optimization of a nonlinear objective function
  with nonlinear equality and/or inequality
  constraints
• Standard NLP form
• System must have more variables than equality
  constraints (under-determined problem)
• Common solvers CONOPT, NPSOL
• Non-convex problems can converge to a local
  optimum

x vector of variables to be determined
(decision variables) h(x) vector function of
equality constraints g(x) vector function of
inequality constraints f(x) scalar objective
function
10
Matlab NLP Solvers lsqnonlin and fmincon

• Nonlinear least-squares lsqnonlin
• x lsqnonlin(_at_fun,x0,lb,ub)
• where fun is a user-defined function that
  returns the vector value F(x),
• x0 is the initial guess (starting point), and
  lb and ub are the bounds on x
• Constrained nonlinear multivariable optimization
  fmincon
• where x, b, beq, lb, and ub are vectors, A
  and Aeq are matrices, c(x)
• and ceq(x) are functions that return vectors,
  and f(x) is a function that
• returns a scalar
• x fmincon(_at_fun,x0,A,b,Aeq,beq,lb,ub,_at_cfun)
• where fun is the function for f(x) and cfun is a
  function that returns c(x) and ceq(x)
• f fun(x) c,ceq cfun(x)

11
Batch Fermentation Example

• Parameter estimation problem for penicillin
  fermentation
• Model equations
• Batch cell growth is modeled by the logistic law
• where y1 is the cell concentration, k1 is the
  growth constant k2 is the cessation (limiting
  nutrient) constant
• Penicillin production is modeled as
• where y2 is the penicillin concentration, k3 is
  the production constant k4 is the degradation
  (hydrolysis) constant
• Dynamic parameter estimation
• Use experimental data from two batch penicillin
  fermentations
• Find values for the unknown parameters (k1, k2,
  k3, k4) that minimize the sum of squared errors
  between the data model predictions

12
Matlab Exercise Batch Data Sets
13
Matlab Exercise Solution

• Load plot the experimental data
• Choose an initial guess, integrate the model,
  plot the simulated profiles
• Estimate parameter values that minimize the sum
  of squared errors between the experimental
  measurements model predictions

gtgt pendat xlsread('penicillin.xls') gtgt tdat
pendat(,1) gtgt ydat pendat(,2end) gtgt
plot(tdat,ydat,'o') gtgt xlabel('Time h') gtgt
ylabel('Concentration')
gtgt k0 0.1 4 0.01 0.01 gtgt y0 0.29 0 gtgt
ts min(tdat) max(tdat) gtgt dy _at_(t,y,k)
k(1)y(1)(1-y(1)/k(2)) k(3)y(1)-k(4)y(2) gtgt
tsim,ysim ode45(dy,ts,y0,,k0) gtgt hold on,
plot(tsim,ysim,'')
gtgt options optimset('Display','iter') gtgt k
lsqnonlin(_at_simerr,k0,,,options,dy,ts,y0,tdat,y
dat) gtgt tsim,ysim ode45(dy,ts,y0,,k) gtgt
plot(tsim,ysim)
14
Matlab Exercise simerr.m
function e simerr(k0,dy,ts,y0,tdat,ydat)
Integrate the model sol ode45(dy,ts,y0,,k0)
Evaluate solution at the data points y
deval(sol,tdat)' Error between data and
model e ydat - y Find missing
measurements n find(isnan(ydat)) Zero error
for missing measurements if isempty(n) e(n)
zeros(size(n)) end