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
