# Optimisation of batch time using genetic algorithm and ODE solver

4 views (last 30 days)
Elan Hutchinson on 14 Mar 2023
Commented: Elan Hutchinson on 15 Mar 2023
Hey guys,
I am trying to adapt the below code to minimise the batch time whilst maintaining the crystal size achieved in the optimisation code below. I am not sure how to go about it given the batch time will change and how to define the time intervals. Any help would be much appreciated.
%Defining temperature bounds
ub = 315*ones(1,5);
lb = 273*ones(1,5);
Aeq = [];
beq =[];
%Defining only cooling allowed (With known initial and final temp)
To = 315;
Tf = 273;
A = [1 0 0 0 0; -1 1 0 0 0; 0 -1 1 0 0; 0 0 -1 1 0; 0 0 0 -1 0];
b = [To 0 0 0 -Tf];
%Other required definitions for GA optimiser
nvars = 5;
fitfun = @Q4;
C0 = 0.0256;
% Opitimiser
options = optimoptions(options,'PlotFcn',{@gaplotbestf}, ...
'Display','iter');
[x,fval] = ga(fitfun,nvars,A,b,Aeq,beq,lb,ub,[],options);
%Initial conditions
x0= [1e-10 0 0 0 C0 0 To];
% Simulation time and time step
t(1) = 0;
tf = 80;
Dt = tf./5;
% Creating matrix to enable plotting of results
Time = [];
XX=[];
%Rate of cooling and discritised time
[f] = @(R) (R>=-2 & R<=0);
for i=1:1:5
t(2) = t(1) + Dt;
x(end,5) = 273;
if i == 1
R = (x(i)-To)./Dt;
@(R) (R>=-2 & R<=0);
else
R = (x(i)-x(i-1))./Dt;
@(R) (R>=-2 & R<=0);
end
% Calling the simulator
[time,X] = ode45(@Q4Model, [t(1) t(2)], x0, [], R);
x0 = X(end,:);
Time =[Time;time];
XX =[XX;X];
t(1) = t(2);
X=real(X);
XX=real(XX);
end
%Plotting Figures
figure
subplot(2,2,1),plot(Time,XX(:,5),'m'), title('a) Concentration'), xlabel('Time (mins)'), ylabel('Concentration (g/g)')
subplot(2,2,2),plot(Time,XX(:,6),'m'), title('b) Mean crystal size'), xlabel('Time (mins)'), ylabel('Mean crystal size (m)')
subplot(2,2,3),plot(Time,XX(:,7),'m'), title('c) Temperature'), xlabel('Time (mins)'), ylabel('Temperature (K)')
subplot(2,2,4),plot(Time,gradient(XX(:,7),Time),'m'), title('d) Rate of cooling'), xlabel('Time (mins)'), ylabel('Cooling rate (K/min)')
function [f] = Q4(T)
% Initial conditions
To = 315;
Tf = 273;
C0 = 0.0256;
x0= [1e-10 0 0 0 C0 0 To];
t(1) = 0;
tf = 80;
% Discretized time
Dt = tf./5;
@(R) (R>=-2 & R<=0);
for i=1:1:5
t(2) = t(1) + Dt;
@(R) (R>=-2 & R<=0);
if i == 1
R = (T(i)-To)./Dt;
else
R = (T(i)-T(i-1))./Dt;
@(R) (R>=-2 & R<=0);
end
% Call the simulator
[t,X] = ode45(@Q4Model, [t(1) t(2)], x0, [], R);
x0 = X(end,:);
t(1) = t(2);
X(end,5)=Tf;
X=real(X);
@(R) (R>=-2 & R<=0);
end
%Objective function
f = -X(end,6);
end
%Mathematical model
function dxdt = Q4Model(t,x,R)
%Variables
mu0 = x(1);
mu1 = x(2);
mu2 = x(3);
mu3 = x(4);
C = x(5);
d = x(6);
T = x(7);
%Constants
kb = 7.8e19;
kg = 0.017;
kv = 0.24;
rho = 1.296e6;
tf=80;
Dt = tf./5;
b=6.2;
g=1.5;
%Concentration Equations
Cs = 1.5846e-5 * T.^2 - 9.0567e-3 * T + 1.3066;
B = kb * (C - Cs).^b;
G = kg * (C - Cs).^g;
%Diffrential Equations
dmu0dt = B;
dmu1dt = G * mu0;
dmu2dt = 2 * G * mu1;
dmu3dt = 3 * G * mu2;
dCdt = -3 * kv * rho * mu3 * G;
dddt = (mu1./mu0)./Dt;
dTdt=R;
dxdt = [dmu0dt; dmu1dt; dmu2dt; dmu3dt; dCdt; dddt; dTdt];
end

Alan Weiss on 14 Mar 2023
I do not know what the batch time and crystal size mean in your model, so I cannot tell you directly what to do.
I can tell you that you are likely wasting a good deal of time in your process by using ga as the optimizer. I see nothing integer-constrained or discontinuous in your model. I thnk that you would do much better to use fmincon as the optimizer, or at the very least patternsearch instead of ga. For fmincon you might need to set a larger-than-default value of the FiniteDifferenceStepSize option; maybe try 1e-6 or 1e-5 to start. See Optimize ODEs in Parallel and Optimizing a Simulation or Ordinary Differential Equation.
Alan Weiss
MATLAB mathematical toolbox documentation
Elan Hutchinson on 15 Mar 2023
I am attempting to minimise the batch time tf in the code below but the solver wont run using ga or fmincon. Think the problem is with defining the time step. If you can see a issue in what I have below that would be super helpful.
%Set up of script
clc
clear
close all
%Defining temperature bounds
ub = 315*ones(1,5);
lb = 273*ones(1,5);
Aeq = [];
beq =[];
%Defining only cooling allowed (With known initial and final temp)
To = 315;
Tf = 273;
A = [1 0 0 0 0; -1 1 0 0 0; 0 -1 1 0 0; 0 0 -1 1 0; 0 0 0 -1 0];
b = [To 0 0 0 -Tf];
%Other required definitions for GA optimiser
nvars = 5;
fitfun = @Q4;
C0 = 0.0256;
% Opitimiser
options = optimoptions(options,'PlotFcn',{@gaplotbestf}, ...
'Display','iter');
[x,fval] = ga(fitfun,nvars,A,b,Aeq,beq,lb,ub,[],options);
%Initial conditions
x0= [1e-10 0 0 0 C0 0 To 80];
% Simulation time and time step
t(1) = 0;
Dt = X(end,8)/5;
% Creating matrix to enable plotting of results
Time = [];
XX=[];
%Rate of cooling and discritised time
[f] = @(R) (R>=-2 & R<=0);
for i=1:1:5
t(2) = t(1) + Dt;
x(end,5) = 273;
if i == 1
R = (x(i)-To)./Dt;
@(R) (R>=-2 & R<=0);
else
R = (x(i)-x(i-1))./Dt;
@(R) (R>=-2 & R<=0);
end
% Calling the simulator
[time,X] = ode45(@Q4Model, [t(1) t(2)], x0, [], R);
x0 = X(end,:);
Time =[Time;time];
XX =[XX;X];
t(1) = t(2);
tf= t(1)+t(2)+t(3)+t(4)+t(5)
X=real(X);
XX=real(XX);
end
%Plotting Figures
figure
subplot(2,2,1),plot(Time,XX(:,5),'m'), title('Concentration'), xlabel('Time (mins)'), ylabel('Concentration (g/g)')
subplot(2,2,2),plot(Time,XX(:,6),'m'), title('Mean crystal size'), xlabel('Time (mins)'), ylabel('Mean crystal size (m)')
subplot(2,2,3),plot(Time,XX(:,7),'m'), title('Temperature'), xlabel('Time (mins)'), ylabel('Temperature (K)')
subplot(2,2,4),plot(Time,XX(:,10),'m'), title('Batch time'), xlabel('Time (mins)'), ylabel('Temperature (K)')
function [f1, f2] = Q4(T)
% Initial conditions
To = 315;
Tf = 273;
C0 = 0.0256;
x0= [1e-10 0 0 0 C0 0 To 1e-10];
t(1) = 0;
% Discretized time
[f] = @(R) (R>=-2 & R<=0);
for i=1:1:5
t(2) = t(1) + Dt;
x(end,5) = 273;
if i == 1
R = (x(i)-To)./Dt;
@(R) (R>=-2 & R<=0);
else
R = (x(i)-x(i-1))./Dt;
@(R) (R>=-2 & R<=0);
end
% Call the simulator
[t,X] = ode45(@Q4Model, [t(1) t(2)], x0, [], R);
x0 = X(end,:);
t(1) = t(2);
X(end,5)=Tf;
X=real(X);
@(R) (R>=-2 & R<=0);
end
%Objective function
f1 = X(end,8);
f2 = -X(end,6)
end
%Mathematical model
function dxdt = Q4Model(t,x,R)
%Variables
mu0 = x(1);
mu1 = x(2);
mu2 = x(3);
mu3 = x(4);
C = x(5);
d = x(6);
T = x(7);
t = x(8);
%Constants
kb = 7.8e19;
kg = 0.017;
kv = 0.24;
rho = 1.296e6;
tf=80;
Dt = tf./5;
b=6.2;
g=1.5;
%Concentration Equations
Cs = 1.5846e-5 * T.^2 - 9.0567e-3 * T + 1.3066;
B = kb * (C - Cs).^b;
G = kg * (C - Cs).^g;
%Diffrential Equations
dmu0dt = B;
dmu1dt = G * mu0;
dmu2dt = 2 * G * mu1;
dmu3dt = 3 * G * mu2;
dCdt = -3 * kv * rho * mu3 * G;
dddt = (mu1./mu0)./Dt;
dTdt=R;
dtfdt= 1
dxdt = [dmu0dt; dmu1dt; dmu2dt; dmu3dt; dCdt; dddt; dTdt; dtfdt];
end

### Categories

Find more on Genetic Algorithm in Help Center and File Exchange

R2022b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!