Scribd
Upload
13 views3 pages

A1 Assignment Part 2

The document outlines a MATLAB code for optimizing a set of parameters in a differential equation model using the lsqnonlin solver. It initializes parameters, defines the ODE system, and sets up an optimization problem to minimize the sum of squared differences between the true solution and the model's output. The optimization results in a set of parameter values that fit the model to the data with a final sum of squares of 0.1686.

Uploaded by

Sreya Banerjee
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
13 views3 pages

A1 Assignment Part 2

The document outlines a MATLAB code for optimizing a set of parameters in a differential equation model using the lsqnonlin solver. It initializes parameters, defines the ODE system, and sets up an optimization problem to minimize the sum of squared differences between the true solution and the model's output. The optimization results in a set of parameter values that fit the model to the data with a final sum of squares of 0.1686.

Uploaded by

Sreya Banerjee
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

Code for Optimisation:

clc
clear all
close all
%initial guess k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14
rtrue = [0.01 0.02 0.05 0.05 0.02 0.02 0.045 0.01 0.02 0.02 0.05 0.45 0.2 0.4];
%first plot points for concentrations at t=0
y0 = [0.492 0.48 0.084 0.101 0.202 0.093];
% time span
tspan = linspace(0,100);
soltrue = ode45(@(t,y)diffun(t,y,rtrue),tspan,y0); %true solution
yvalstrue = deval(soltrue,tspan); % solution vectors returned
%plot
for i = 1:6
subplot(3,2,i)
plot(tspan,yvalstrue(i,:))
title(['y(',num2str(i),')'])
end
% framing optimisation
r = optimvar('r',14,"LowerBound",0.001,"UpperBound",10); %14 r values where
0.01<=r<=10
myfcn = fcn2optimexpr(@RtoODE,r,tspan,y0); % optimisation function RtoODE
obj = sum(sum((myfcn - yvalstrue).^2)); %objective function as the sum of squared
differences between the ODE solution and the solution with true parameters
prob = optimproblem("Objective",obj);
show(prob)
r0.r = [1 1 1 1 1 1 1 1 1 1 1 1 1 1]; % initial guess to call solver
[rsol,sumsq] = solve(prob,r0)
disp(rsol.r)

function solpts = RtoODE(r,tspan,y0)


sol = ode45(@(t,y)diffun(t,y,r),tspan,y0);
solpts = deval(sol,tspan);
end

function dydt = diffun(~,y,r)


dydt = zeros(6,1);
dydt(1) = -r(1)*y(1)+r(2)*y(3);
dydt(2) = -r(3)*y(2)+r(4)*y(3)-r(11)*y(2)+r(12)*y(5)-r(13)*y(2)+r(14)*y(6);
dydt(3) = -r(5)*y(3) + r(6)*y(4) - r(2)*y(3) +r(1)*y(1) +r(4)*y(2) - r(3)*y(3);
dydt(4) = r(8)*y(5)-r(7)*y(4)+r(5)*y(3)-r(6)*y(4);
dydt(5) = r(7)*y(4)-r(8)*y(5)+r(10)*y(6)-r(9)*y(5)+r(11)*y(2)-r(12)*y(5);
dydt(6) = r(9)*y(5)-r(10)*y(6)+r(13)*y(2)-r(14)*y(6);
end
Output:

OptimizationProblem :

Solve for:
r

minimize :
sum(sum((RtoODE(r, extraParams{1}, extraParams{2}) -
extraParams{3}).^2, 1))

extraParams

variable bounds:
0.001 <= r(1) <= 10
0.001 <= r(2) <= 10
0.001 <= r(3) <= 10
0.001 <= r(4) <= 10
0.001 <= r(5) <= 10
0.001 <= r(6) <= 10
0.001 <= r(7) <= 10
0.001 <= r(8) <= 10
0.001 <= r(9) <= 10
0.001 <= r(10) <= 10
0.001 <= r(11) <= 10
0.001 <= r(12) <= 10
0.001 <= r(13) <= 10
0.001 <= r(14) <= 10

Solving problem using lsqnonlin.

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to
its initial value is less than the value of the function tolerance.

<stopping criteria details>


rsol =

struct with fields:

r: [14×1 double]

sumsq =

0.1686

1.1296
2.1450
0.2951
0.2952
0.6858
1.6413
0.4987
0.2491
7.8553
1.1879
0.5648
0.0010
1.0630
3.5602

You might also like