why didn't the event function detect the events?

조회 수: 5 (최근 30일)
汉武 沈
汉武 沈 2021년 1월 21일
댓글: 汉武 沈 2021년 1월 21일
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

답변 (2개)

Mischa Kim
Mischa Kim 2021년 1월 21일
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
汉武 沈
汉武 沈 2021년 1월 21일
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!

댓글을 달려면 로그인하십시오.


Jan
Jan 2021년 1월 21일
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
汉武 沈
汉武 沈 2021년 1월 21일
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!

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

태그

제품

Community Treasure Hunt

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

Start Hunting!

Translated by