possibility to add a timeror on ODE15s which after a certain time limit, the solver is stopped or gives a null output

8 views (last 30 days)
Shreenath Krishnamurthy
Shreenath Krishnamurthy on 18 Mar 2020
Answered: Shreenath Krishnamurthy on 6 Apr 2020
Hi I am using genetic algorithm gamultiobj to minimize an objective function (given in mainfile). The objective function (given in fobj) consists of differential algebraic equations. For certain set of the parameters, the odes take day or two to solve, which is a long time.
Is there way to limit ode15s max to only take about 15 minutes , if it exceeds 15 minutes it should stop and take a new input from the GA. The functions i use are given below
mainfile
%% Main file for MOO of 6 step cycle. Obj 1: Max Pu, Obj 2: Max En
% clc
% clear
global iteration
iteration=0;
t1 = tic;
NumVar=7;
Xlow=[10 1 10 1 0.1 0.1 0.1];
Xup =[200 100 200 100 0.5 0.3 2.0];
A=[0 0 0 0 -1 1 0];
b=[-0.02];
options = optimoptions('gamultiobj','UseParallel', true,'UseVectorized', false,'ParetoFraction',0.35,...
'PopulationSize',140,'Generations',30, 'CrossoverFraction',0.5,'CrossoverFcn',@crossoverintermediate ,'MutationFcn',...
@mutationadaptfeasible,'InitialPopulation', X0);%'outputfcns',@gaoutput,...
parpool(8);
[X,fval,exitflag,output,final_pop]=gamultiobj(@fobj,NumVar,A,b,[],[],...
Xlow,Xup,options);
function [f]=fobj(X)
global iteration
iteration=0.0;
iteration=iteration+1;
persistent P1 P2 P3 P4
%parameters
param1=X(1);
param2=X(2);
param3=X(3);
param4=X(4);
param5=X(5);
param6=X(6);
param7=X(7);
%This sections specifies the constanta
%==============================================================================================
%Spatial dimensions and initial grid
ngrid=30;
h=1/ngrid;
ic1=ones(1,ngrid);
ic2=ones(1,ngrid);
ic3=ones(1,ngrid);
ic4=zeros(1,ngrid);
ic5=ones(1,ngrid);
ic6=zeros(1,ngrid);
ic=[ic1 ic2 ic3 ic4 ic5 ic6];
Y=ic;
%save RESULTS_shre1
try
while swtch==0
n=n+1;
if n<=nswitch
%%step 1
Y=YALL(end,1:6*ngrid);
Y(6*ngrid+1)=PL/PH;Y(6*ngrid+2)=0;Y(6*ngrid+3)=yaf; Y(6*ngrid+4)=1;Y(6*ngrid+5)=10;Y(6*ngrid+6)=10*yaf;
vinit=0;
delt=0.01;
options=odeset('RelTol', [1e-3], 'AbsTol', [1e-4]);
t1=cputime;
[time,YALL]=odeq15s(@f1, 0:delt:t1,Y,options);
YALL1=YALL;
%% step 2
Y=YALL(end,1:6*ngrid);
Y(6*ngrid+(1:2))=0;
vinit=0;
delt=0.1;
options=odeset('RelTol', [1e-3], 'AbsTol', [1e-3],'Vectorized','on');
t1=cputime;
[time,YALL]=ode15s(@f2, 0:delt:t2,Y,options);
YALL2=YALL;
%% step 3
Y=YALL(end,1:6*ngrid);
Y(6*ngrid+(1:3))=0;
vinit=0;
delt=0.1;
options=odeset('RelTol', [1e-4], 'AbsTol', [1e-4],'Vectorized','on');
t1=cputime;
[time,YALL]=ode15s(@f3, 0:delt:t3,Y,options);
YALL3=YALL;
%% step4
Y=YALL(end,1:6*ngrid);
for i=6:-1:1
abc((i-1)*ngrid+(1:ngrid))=Y(i*ngrid:-1:(i-1)*ngrid+1);
end
clear Y; Y=abc;clear abc;
vinit=0;
delt=0.1;
options=odeset('RelTol', [1e-4], 'AbsTol', [1e-4],'Vectorized','on');
t1=cputime;
[time,YALL]=ode15s(@f4, 0:delt:t4,Y,options);
for i=6:-1:1
YALL(:,(i-1)*ngrid+(1:ngrid))=YALL(:,i*ngrid:-1:(i-1)*ngrid+1);
end
YALL4=YALL;
%%step5
Y=YALL(end,1:6*ngrid);
for i=6:-1:1
abc((i-1)*ngrid+(1:ngrid))=Y(i*ngrid:-1:(i-1)*ngrid+1);
end
clear Y; Y=abc;clear abc;
Y(6*ngrid+(1:4))=0;
vinit=0;
delt=0.1;
options=odeset('RelTol', [1e-4], 'AbsTol', [1e-4],'Vectorized','on');
t1=cputime;
[time2,YALL]=ode15s(@f5, 0:delt:t5,Y,options);
for i=6:-1:1
YALL(:,(i-1)*ngrid+(1:ngrid))=YALL(:,i*ngrid:-1:(i-1)*ngrid+1);
end
YALL7=YALL;
end
if n>=nswitch
%%step 6
Y=YALL(end,1:6*ngrid);
Y(6*ngrid+(1:6))=0;
for i=6:-1:1
abc((i-1)*ngrid+(1:ngrid))=Y(i*ngrid:-1:(i-1)*ngrid+1);
end
clear Y; Y=abc;clear abc;
[time,YALL]=ode15s(@f6, 0:delt:t6,Y,options);
t2press=cputime-t1;
for i=6:-1:1
YALL(:,(i-1)*ngrid+(1:ngrid))=YALL(:,i*ngrid:-1:(i-1)*ngrid+1);
end
YALL5=YALL;
fluegasflag=0.0;
if Pstar<0.99
fluegasflag=1.0;
Y=YALL(end,1:6*ngrid);
Y(6*ngrid+(1:6))=0;
delt=0.01;
options=odeset('RelTol', [1e-3], 'AbsTol', [1e-4],'Vectorized','on');
t1=cputime;
[time,YALL]=ode15s(@f7, 0:delt:t7,Y,options);YALL6=YALL;
t2press=cputime-t1;
%% step 2
Y=YALL(end,1:6*ngrid);
Y(6*ngrid+(1:2))=0;
vinit=0;
delt=0.1;
options=odeset('RelTol', [1e-3], 'AbsTol', [1e-3],'Vectorized','on');
t1=cputime;
[time,YALL]=ode15s(@f2, 0:delt:t2,Y,options);
YALL2=YALL;
%% step 8
Y=YALL(end,1:6*ngrid);
Y(6*ngrid+(1:2))=0;
vinit=0;
delt=0.1;
options=odeset('RelTol', [1e-3], 'AbsTol', [1e-3],'Vectorized','on');
t1=cputime;
[time,YALL]=ode15s(@f8, 0:delt:t2,Y,options);
YALL2=YALL;
%% step 3
Y=YALL(end,1:6*ngrid);
Y(6*ngrid+(1:3))=0;
vinit=0;
delt=0.1;
options=odeset('RelTol', [1e-4], 'AbsTol', [1e-4],'Vectorized','on');
t1=cputime;
[time,YALL]=ode15s(@f3, 0:delt:t3,Y,options);
YALL3=YALL;
%% step4
Y=YALL(end,1:6*ngrid);
for i=6:-1:1
abc((i-1)*ngrid+(1:ngrid))=Y(i*ngrid:-1:(i-1)*ngrid+1);
end
clear Y; Y=abc;clear abc;
vinit=0;
delt=0.1;
options=odeset('RelTol', [1e-4], 'AbsTol', [1e-4],'Vectorized','on');
t1=cputime;
[time,YALL]=ode15s(@f4, 0:delt:t4,Y,options);
for i=6:-1:1
YALL(:,(i-1)*ngrid+(1:ngrid))=YALL(:,i*ngrid:-1:(i-1)*ngrid+1);
end
YALL4=YALL;
%%step5
Y=YALL(end,1:6*ngrid);
for i=6:-1:1
abc((i-1)*ngrid+(1:ngrid))=Y(i*ngrid:-1:(i-1)*ngrid+1);
end
clear Y; Y=abc;clear abc;
Y(6*ngrid+(1:4))=0;
vinit=0;
delt=0.1;
options=odeset('RelTol', [1e-4], 'AbsTol', [1e-4],'Vectorized','on');
t1=cputime;
[time2,YALL]=ode15s(@f5, 0:delt:t5,Y,options);
for i=6:-1:1
YALL(:,(i-1)*ngrid+(1:ngrid))=YALL(:,i*ngrid:-1:(i-1)*ngrid+1);
end
YALL7=YALL;
end
%% Nett calculations
if n<=nswitch
P1=....
P2= ....
P3=....
P4= ....
end
if n>nswitch
P1=....
P2= ....
P3=....
P4= ....
end
end
%% Cyclic steady state condition
if n==300
swtch=1;
end
if n>50 % will check after 50 cycles
ind=find(masslist(end-5+1:end)<0.5);
if length(ind)==5
swtch=1;
end
end
end
TolPur=0.95;
TolRec=0.9;
f=zeros(2,1);
f(1)= F1(P1,P2,P3,P4)
f(2)= F1(P1,P2,P3,P4);
catch
f = [200000;200000];
end
filename = [tempname('DIRECTORY'),'.txt'];
%tmpName = [filename,'.txt'];
fid = fopen(filename,'w');
fprintf(fid, '%f %f %f %f %f %f %f %f %f %f %f\n',X(1),X(2),X(3),X(4),X(5),X(6),X(7),purity,recovery,Prod,Energy);
fclose(fid);
%==================================================================================================================
function dydt=f(time,y)
dydt=zeros(6*ngrid+6,size(y,2));
dydt(6*ngrid+1,:)=(1-PL/PH)*(ATP*L/V0)*exp(-ATP*time*L/V0);
dydt (1:ngrid,:) =.....
dydt(ngrid+1:2*ngrid,:)=....
dydt(2*ngrid+1:3*ngrid,:)= ....
dydt (3*ngrid+1:4*ngrid,:) =.....
dydt(4*ngrid+1:5*ngrid,:)=....
dydt(5*ngrid+1:6*ngrid,:)= ....
;
end
function dydt=f2(time,y)
dydt=zeros(6*ngrid+2,size(y,2));
ve(ngrid,size(y,2))=0;
ve0=1;
h=1/ngrid;
dydt (1:ngrid,:) =.....
dydt(ngrid+1:2*ngrid,:)=....
dydt(2*ngrid+1:3*ngrid,:)= ....
dydt (3*ngrid+1:4*ngrid,:) =.....
dydt(4*ngrid+1:5*ngrid,:)=....
dydt(5*ngrid+1:6*ngrid,:)= ....
;
end
function dydt=f3(time,y)
dydt=zeros(6*ngrid+3,size(y,2));
dydt (1:ngrid,:) =.....
dydt(ngrid+1:2*ngrid,:)=....
dydt(2*ngrid+1:3*ngrid,:)= ....
dydt (3*ngrid+1:4*ngrid,:) =.....
dydt(4*ngrid+1:5*ngrid,:)=....
dydt(5*ngrid+1:6*ngrid,:)= ....
;
end
function dydt=f4(time,y)
dydt=zeros(6*ngrid+3,size(y,2));
dydt (1:ngrid,:) =.....
dydt(ngrid+1:2*ngrid,:)=....
dydt(2*ngrid+1:3*ngrid,:)= ....
dydt (3*ngrid+1:4*ngrid,:) =.....
dydt(4*ngrid+1:5*ngrid,:)=....
dydt(5*ngrid+1:6*ngrid,:)= ....
;
end
function dydt=f5(time,y)
dydt=zeros(6*ngrid+6,size(y,2));
dydt(6*ngrid+1,:)=(1-Pstar)*(ATP*L/V0)*exp(-ATP*time*L/V0);
dydt (1:ngrid,:) =.....
dydt(ngrid+1:2*ngrid,:)=....
dydt(2*ngrid+1:3*ngrid,:)= ....
dydt (3*ngrid+1:4*ngrid,:) =.....
dydt(4*ngrid+1:5*ngrid,:)=....
dydt(5*ngrid+1:6*ngrid,:)= ....
;
end
function dydt=fLPP(time,y)
dydt=zeros(6*ngrid+6,size(y,2));
dydt (1:ngrid,:) =.....
dydt(ngrid+1:2*ngrid,:)=....
dydt(2*ngrid+1:3*ngrid,:)= ....
dydt (3*ngrid+1:4*ngrid,:) =.....
dydt(4*ngrid+1:5*ngrid,:)=....
dydt(5*ngrid+1:6*ngrid,:)= ....
;
end
function dydt=f7(time,y)
%=============================================================================================================
dydt=zeros(6*ngrid+4,size(y,2));
dydt (1:ngrid,:) =.....
dydt(ngrid+1:2*ngrid,:)=....
dydt(2*ngrid+1:3*ngrid,:)= ....
dydt (3*ngrid+1:4*ngrid,:) =.....
dydt(4*ngrid+1:5*ngrid,:)=....
dydt(5*ngrid+1:6*ngrid,:)= ....
;
end
function dydt=f8(time,y)
%=============================================================================================================
dydt=zeros(6*ngrid+4,size(y,2));
dydt (1:ngrid,:) =.....
dydt(ngrid+1:2*ngrid,:)=....
dydt(2*ngrid+1:3*ngrid,:)= ....
dydt (3*ngrid+1:4*ngrid,:) =.....
dydt(4*ngrid+1:5*ngrid,:)=....
dydt(5*ngrid+1:6*ngrid,:)= ....
;
end
function [value,isterminal,direction] = events(time,y)
value=y(5*ngrid+1)-0.99;
isterminal = 1;
direction=1;
end
end
time = toc(t1);
delete(gcp('nocreate')) %to close parpool

Sign in to answer this question.

Answers (2)

Shreenath Krishnamurthy
Shreenath Krishnamurthy on 6 Apr 2020
@Josh Meyer. Thanks for the reply. I also found this link which is very similar
https://se.mathworks.com/matlabcentral/answers/36683-timeout-using-parfor-loop-and-ode15s
I have tried both and they seem to work. I agree with you that this is not the most efficient way unlike Jpattern, but it is some how doing the job.
Josh Meyer
Josh Meyer on 3 Apr 2020
Edited: Josh Meyer on 3 Apr 2020
You could implement an event function to do that:
https://www.mathworks.com/help/matlab/math/ode-event-location.html
Use startTime = datetime('now') to find the current time in your script, then pass that value to the event function and add an if-statement that triggers a terminal event when the timer gets to the point you want:
[value,isterminal,direction] = myTimer(t,y,startTime)
stopTime = startTime + minutes(15);
isterminal = 1;
direction = 0;
if datetime('now') >= stopTime
value = 0;
else
value = 1;
end
end
One upside to this is that when the event triggers, ode15s returns the solution and current time, so you can adjust your inputs and then restart the integration from that point if you want (the ballode example shows stopping/restarting with an event function).
That said, the downside to this is that each step in the solution checks what time it is! You'd be able to terminate after a specified amount of time has passed, but the time spent solving will be less efficient.

Categories

Find more on Function Creation in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!