my code doesn't work as it should,need help

Hello, I am working on a project and have run into some issues,namely I am trying to make a certain model and to do so first I need to solve system of 3 differential equations.Inputs need to change over time to simulate changing real life conditions,so I made matrices with some realistic values and for loop.When looking at graph ,I realised that only last values in input matrices alter values of temperatures in graph. So here is the code.

Wc = [1 ;2 ;3 ;0 ;0 ;1 ;2 ;0 ;0 ;1]; 
Fi_s = [0; 30; 70; 0; 0; 20; 90; 20; 0; 15]; 
Ta = [0; 3; 7; 5; 2; 4; 6; 3; 0; 2];
inputs = [Wc Fi_s Ta]
T_0 = [10 10 10];
T_sim = T_0;
Ts = 48*3600; 
x = zeros(3,1);
%x should be vector with x(1) being Tr, x(2) being Tf %and x(3) being Tw %so x = [Tr;Tf;Tw]
for k=1:10
    [t, x] = ode23(@(t, x) model_of_pump(t, x, inputs(k,:)), [0 Ts], T_0); 
    T_0=x(end,:); 
    T_sim=[T_sim;T_0];
end
figure 
plot(t/3600,x(:,1),'-r',t/3600,x(:,2),'-b',t/3600,x(:,3),'-g') legend('Tr', 'Tf', 'Tw') 
xlabel('t') 
ylabel('Temperature')

and the function code

function dxdt = model_of_pump(t,x,inputs)
p = 0.1; ni = 3;
Wc = inputs(1); Fi_s = inputs(2); Ta = inputs(3);
dxdt = zeros(3,1);
Qra = 28 * (x(1,end) - Ta);
Qfr = 624 * (x(2,end) - x(1,end)); 
Qwf = 28 * (x(3,end) - x(2,end)); Qc = ni * Wc;
dxdt(1) = (Qfr - Qra + (1 - p) * Fi_s)/810;
dxdt(2) = (Qwf - Qfr + p * Fi_s)/3315;
dxdt(3) = (Qc - Qwf)/836;
end

[Kevin Holly]

How did you want to display all 10 inputs? I placed a hold on below to graph each iteration of k.

Wc = [1 ;2 ;3 ;0 ;0 ;1 ;2 ;0 ;0 ;1]; 
Fi_s = [0; 30; 70; 0; 0; 20; 90; 20; 0; 15]; 
Ta = [0; 3; 7; 5; 2; 4; 6; 3; 0; 2];
inputs = [Wc Fi_s Ta]

inputs = 10×3

1 0 0 2 30 3 3 70 7 0 0 5 0 0 2 1 20 4 2 90 6 0 20 3 0 0 0 1 15 2

T_0 = [10 10 10];
T_sim = T_0;
Ts = 48*3600; 
x = zeros(3,1);
%x should be vector with x(1) being Tr, x(2) being Tf %and x(3) being Tw %so x = [Tr;Tf;Tw]

figure 
hold on
for k=1:10
    [t, x] = ode23(@(t, x) model_of_pump(t, x, inputs(k,:)),[0 Ts], T_0); 
    T_0=x(end,:); 
    T_sim=[T_sim;T_0];
    plot(t/3600,x(:,1),'-r',t/3600,x(:,2),'-b',t/3600,x(:,3),'-g')
end

legend('Tr', 'Tf', 'Tw') 
xlabel('t') 
ylabel('Temperature')

function dxdt = model_of_pump(t,x,inputs)
    p = 0.1; ni = 3;
    Wc = inputs(1); Fi_s = inputs(2); Ta = inputs(3);
    dxdt = zeros(3,1);
    Qra = 28 * (x(1,end) - Ta);
    Qfr = 624 * (x(2,end) - x(1,end)); 
    Qwf = 28 * (x(3,end) - x(2,end)); Qc = ni * Wc;
    dxdt(1) = (Qfr - Qra + (1 - p) * Fi_s)/810;
    dxdt(2) = (Qwf - Qfr + p * Fi_s)/3315;
    dxdt(3) = (Qc - Qwf)/836;
end

[Lovre]

First thank you for taking your time and answering my question.

What I wanted to do here is simulate work of a heat pump ,with simulation time lasting 2 days , and during that period as conditions change my inputs also change, so for example if my simulation time is 2 days or 48 hours then since my matrix has 10 values I want the first input values to account for first 1/10 of total simulation time or 4.8 hours, and for calculating temperatures for next 4.8 hours I want to use input values with index 2 and so on.What I am trying to get in the end is one graph which would show how change of input arguments changes outputs,so for example there wouldnt be 10 green lines as in your graph but only one green line whose value would change over time.

What worries me the most here is that I noticed how if I change only last value of matrix Wc outputs change drastically and if I change any other value in that matrix or all of them ,outputs stay the same.So I am doubting if I wrote my code in correct way.

To give you an example,if I take this matrix [1 ;2 ;3 ;0 ;0 ;1 ;2 ;0 ;0 ;1] for Wc I get the same outputs as with this matrix [0 ;0 ;0 ;0 ;0 ;0 ;0 ;0 ;0 ;1].And that should not happen.Wc is power supplied to the pump so when value is higher temperature should rise and when there is 0 temperature should slowly decrease.

There is also a question of why temperature increases by such a small amount,because for input power I used room temperature(red line) should be around 20 degrees Celsius.Maybe I tried to solve diff equations in a wrong way(I checked multilpe times if I wrote equations correctly and if I took right constants), but if thats not the case then I should try to investigate why is that so some other time,now my priority is figuring out what's wrong with Wc.

[Walter Roberson]

"and for calculating temperatures for next 4.8 hours I want to use input values with index 2 and so on."

You need to stop the ode45 computation every time there is a discontinuity in the control inputs, and then resume with the boundaries adjusted appropriately. The mathematics of the RK process is only valid if the inputs have continuous second derivatives for the duration of the ode45 call.

Answer 1

[Sam Chak]

Hi @Lovre

Not sure if this is what you want. but I just modified a little to simulate the heat pump in accordance to your requirement (4.8 hours for each input condition)

Ts    = 4.8*3600;   % duration for each interval

Added odeset() to improve the results, and the integration time is changed to from [0 Ts] to [Ts*(k-1) Ts*k].

for k = 1:10
    opt    = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
    [t, x] = ode23(@(t, x) model_of_pump(t, x, inputs(k,:)),[Ts*(k-1) Ts*k], T_0, opt); 
    T_0    = x(end, :); 
    % T_sim  = [T_sim; T_0];
    % plot(t/3600, x(:,1), '-r', t/3600, x(:,2), '-b', t/3600, x(:,3), '-g')
    plot(t/3600, x(:,1), 'linewidth', 1.5, 'Color', [0.4940, 0.1840, 0.5560])
    plot(t/3600, x(:,2), 'linewidth', 1.5, 'Color', [0.4660, 0.6740, 0.1880])
    plot(t/3600, x(:,3), 'linewidth', 1.5, 'Color', [0.9290, 0.6940, 0.1250])
end

[Lovre]

Results I get with your code are close to what I wanted but temperature shouldn't jump like that, namely when for one of the Wc inputs I take 0, temperature shouldn't just drop to zero but slowly decrease so that I could have let's say 4th,5th and 6th input set to 0 and if I put some greater value in 3rd input then my temperatures won't drop to zero immediately on 4th but maybe on 7th iteration step, but then I put some greater value there and temperature rises again.Now temperature shouldn't drop all the way to 0, but should instead drop just a few degrees and go up again.I'l try to adjust inputs a bit to get desired results as it's probably possible to get results I want with this code.

Thank you for help Sam

[Sam Chak]

Hi @Lovre

Technicality-wise, MATLAB is able to produce the multi-interval simulation as requested. On the accuracy or expectation of the results, I may not be able to help because I'm unfamiliar with the heat pump. So, I cannot check whether the differential equations are correctly or not. Like you have suspected. since the physics does not explain the temperature drop, then the culprit might one or more of the following:

• parameters
• signs of terms
• structure of equation