MATLAB Answers

why didn't the event function detect the events?

2 views (last 30 days)
汉武 沈
汉武 沈 on 21 Jan 2021
Commented: 汉武 沈 on 21 Jan 2021
Hi all! I'm solving the thermostat model, which presents the characteristic of hybrid dynamical systems, when the ith (here i=1,2) room's temperature reduces to T_ref-0.5, the thermal is off, when the ith (here i=1,2) room's temperature raises to T_ref+0.5, the thermal is on, I use the ode suits with events, but it failed, plus the program also runs a long time, about 40 seconds(performance of computer is not bad), can anyone help me to fix the problem?
This is ode45 codes, but failed.
tic
%options=odeset('Events',@Events1,'Events',@Events2,'AbsTol',1e-8,'RelTol',1e-8);
options=odeset('AbsTol',1e-8,'RelTol',1e-8);
y0 = [21;25;17];
tspan=0:1000;
[tout,yout]=ode45(@Tq_Tj,tspan, y0,options);
figure(1)
plot(tout,yout(:,1),'k');hold on
plot(tout,yout(:,2),'-.');
plot(tout,yout(:,3),'--');
legend('room1','room2','wall')
toc
function f=Tq_Tj(~,y)
C1=550;C2=600;Cw=580;R1=1.8;R2=2.2;R1w=2;R2w=2;Ta=32;p1=14;p2=14;T_ref=20;
if(y(1)>T_ref+0.5)
p1=14;
elseif(y(1)<T_ref-0.5)
p1=0;
end
if(y(2)>T_ref+0.5)
p2=14;
elseif(y(2)<T_ref-0.5)
p2=0;
end
f=[((Ta-y(1))/R1+(y(3)-y(1))/R1w-p1)/C1; ...
((Ta-y(2))/R2+(y(3)-y(2))/R2w-p2)/C2;...
((y(1)-y(3))/R1w+(y(2)-y(3))/R2w)/Cw];
end
function [g,isterminal,direction]=Events1(~,y)
T_ref=20;
g=[y(1)-(T_ref-0.5);y(1)-(T_ref+0.5)];
isterminal=[0;0];
direction=[0;0];
end
function [g,isterminal,direction]=Events2(~,y)
T_ref=20;
g=[y(2)-(T_ref-0.5);y(2)-(T_ref+0.5)];
isterminal=[0;0];
direction=[0;0];
end
this is another way to solve the problem, you can run it to get the right results.
clc
clear
tic
%为各变量赋初值
delta_q=1e-4;
q1=21;q2=25;q3=17;
x1=q1;x2=q2;x3=q3;
t=0;delta_t=0;
A=zeros(10,6);
n=0;
C1=550;C2=600;Cw=580;R1=1.8;R2=2.2;R1w=2;R2w=2;Ta=32;p1=14;p2=14;m1=1;m2=1;T_ref=20;
%C1=0.56;C2=0.28;Cw=0.2;R1=5;R2=6;R1w=3;R2w=3;Ta=20;p1=4;p2=4;m1=1;m2=1;T_ref=20;
%开始while循环
while (t<1000)
if(x1>T_ref+0.5)
m1=1;
elseif(x1<T_ref-0.5)
m1=0;
end
if(x2>T_ref+0.5)
m2=1;
elseif(x2<T_ref-0.5)
m2=0;
end
Dx1=((Ta-q1)/R1+(q3-q1)/R1w-m1*p1)/C1;
Dx2=((Ta-q2)/R2+(q3-q2)/R2w-m2*p2)/C2;
Dx3=((q1-q3)/R1w+(q2-q3)/R2w)/Cw;
DDx1=(-Dx1/R1+(Dx3-Dx1)/R1w)/C1;
DDx2=(-Dx2/R2+(Dx3-Dx2)/R2w)/C2;
DDx3=((Dx1-Dx3)/R1w+(Dx2-Dx3)/R2w);
DDDx1=(-DDx1/R1+(DDx3-DDx1)/R1w)/C1;
DDDx2=(-DDx2/R2+(DDx3-DDx2)/R2w)/C2;
DDDx3=((DDx1-DDx3)/R1w+(DDx2-DDx3)/R2w);
%求Δt1和Δt2的值
delta_t1=sqrt(2*delta_q/abs(DDx1));
delta_t2=sqrt(2*delta_q/abs(DDx2));
delta_t3=sqrt(2*delta_q/abs(DDx3));
delta_tmin=min([delta_t1 delta_t2 delta_t3]);
%比较Δt1和Δt2的大小,进而确定q1和q2谁先跃迁
if (delta_t1==delta_tmin)
delta_t=delta_t1;
t=t+delta_t;
caribe_1=-0.125*DDDx1*(exp(-2*delta_t)-1)+0.5*(DDx1+0.5*DDDx1)*delta_t^2+(Dx1-0.25*DDDx1)*delta_t;
caribe_2=-0.125*DDDx2*(exp(-2*delta_t)-1)+0.5*(DDx2+0.5*DDDx2)*delta_t^2+(Dx2-0.25*DDDx2)*delta_t;
caribe_3=-0.125*DDDx3*(exp(-2*delta_t)-1)+0.5*(DDx3+0.5*DDDx3)*delta_t^2+(Dx3-0.25*DDDx3)*delta_t;
x1=x1+caribe_1;
x2=x2+caribe_2;
x3=x3+caribe_3;
q1=x1;
q2=q2+caribe_2;
q3=q3+caribe_3;
elseif (delta_t2==delta_tmin)
delta_t=delta_t2;
t=t+delta_t;
caribe_1=-0.125*DDDx1*(exp(-2*delta_t)-1)+0.5*(DDx1+0.5*DDDx1)*delta_t^2+(Dx1-0.25*DDDx1)*delta_t;
caribe_2=-0.125*DDDx2*(exp(-2*delta_t)-1)+0.5*(DDx2+0.5*DDDx2)*delta_t^2+(Dx2-0.25*DDDx2)*delta_t;
caribe_3=-0.125*DDDx3*(exp(-2*delta_t)-1)+0.5*(DDx3+0.5*DDDx3)*delta_t^2+(Dx3-0.25*DDDx3)*delta_t;
x1=x1+caribe_1;
x2=x2+caribe_2;
x3=x3+caribe_3;
q2=x2;
q1=q1+caribe_1;
q3=q3+caribe_3;
elseif (delta_t3==delta_tmin)
delta_t=delta_t3;
t=t+delta_t;
caribe_1=-0.125*DDDx1*(exp(-2*delta_t)-1)+0.5*(DDx1+0.5*DDDx1)*delta_t^2+(Dx1-0.25*DDDx1)*delta_t;
caribe_2=-0.125*DDDx2*(exp(-2*delta_t)-1)+0.5*(DDx2+0.5*DDDx2)*delta_t^2+(Dx2-0.25*DDDx2)*delta_t;
caribe_3=-0.125*DDDx3*(exp(-2*delta_t)-1)+0.5*(DDx3+0.5*DDDx3)*delta_t^2+(Dx3-0.25*DDDx3)*delta_t;
x1=x1+caribe_1;
x2=x2+caribe_2;
x3=x3+caribe_3;
q3=x3;
q1=q1+caribe_1;
q2=q2+caribe_2;
end
n=n+1;
A(n,1)=t;
A(n,2)=x1;
A(n,3)=x2;
A(n,4)=x3;
A(n,5)=m1;
A(n,6)=m2;
end
% A;
figure(2);
plot(A(:,1),A(:,2),'--');hold on
plot(A(:,1),A(:,3),'-.');
plot(A(:,1),A(:,4),'-');
grid on
% plot(A(:,1),A(:,5)*22);hold on
% plot(A(:,1),A(:,6)*22);
legend('room1','room2','wall')
toc

Answers (2)

Mischa Kim
Mischa Kim on 21 Jan 2021
In your ode45 call you are not returning any event information. Replace with
[tout,yout,te,ye,ie] = ode45(@Tq_Tj,tspan, y0,options);
and you should receive events information such as event times te, solution value at each of the event times, etc.
Regarding execution time, this also depends on the tolerances you set in the ode options. I do not know what your requirements are, however, try lowering the tolerances (or removing them entirely, just to see what happens) to see if the results are still satisfactory.
  1 Comment
汉武 沈
汉武 沈 on 21 Jan 2021
Hello, I just did what you said, but the results is not satisfactory, the result is below.
The model is in fact that two air conditionings placed in two rooms, the wall is in the medium, I write it down in a paper, hope you can understand it.
And another thing What shoud I do with the 'te' and 'ie' information? Is the event function right, how to correct it?
Thank you very much!

Sign in to comment.


Jan
Jan on 21 Jan 2021
Just my standard interjection: Matlab's ODE integrators are designed to integrate smooth functions. Your function Tq_Tj() conatins discontinuities. This can reduce the step size until the processing time gets huge and the result is dominated by rounding errors.
In
options=odeset('Events',@Events1,'Events',@Events2, ...
you overwrite the Event function. Only the 2nd one is considered. Instead of trying to provide 2 functions, implement both events in 1 function. You have to stop the integration, when an event occurs to avoid the discontinuity. Restart the integration with the final values of the former interval.
  1 Comment
汉武 沈
汉武 沈 on 21 Jan 2021
Hi,though I have understood your meaning, but sadly I don't know how to code, forgive my poor ability :), if you are free, can you help me get it done?
Thank you, anyway!

Sign in to comment.

Tags

Products

Community Treasure Hunt

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

Start Hunting!

Translated by