Optimisation of batch time using genetic algorithm and ODE solver

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(@ga,'MutationFcn',@mutationadaptfeasible);

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

Answer 1

[Alan Weiss]

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]

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(@ga,'MutationFcn',@mutationadaptfeasible);

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