why didn't the event function detect the events?

Question

汉武沈 2021-1-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/722808-why-didn-t-the-event-function-detect-the-events

评论：汉武沈 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
       

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mischa Kim 2021-1-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/722808-why-didn-t-the-event-function-detect-the-events#answer_602878

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -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!

请先登录，再进行评论。

Answer 2

Jan 2021-1-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/722808-why-didn-t-the-event-function-detect-the-events#answer_602948

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -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!

请先登录，再进行评论。

why didn't the event function detect the events?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

why didn't the event function detect the events?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论