why didn't the event function detect the events?
1 回表示 (過去 30 日間)
古いコメントを表示
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
0 件のコメント
回答 (2 件)
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.
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.
参考
カテゴリ
Help Center および File Exchange で Ordinary Differential Equations についてさらに検索
製品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!