Solving differential equation in simulink

Hello Everyone

I am solving a differential equation (attached) of greenhouse energy balance on SIMULINK which has basically two unknowns (Tin and qg). One of them (qg) is model output and other (Tin) is also not know from any means. Now I have to calculate both of these parameters from this differential equation. The only thing that I know is the 6 factors on rightside of the equation but all of them are dependent on Tin and there is a operating set point for temperature is given as 20 for day and 16 for night for greenhouse. Can someone help me to solve this equation?

My simulink model is also attached.Plzzzz guide me

[hosein Javan]

if you are solving an ode then you have only one differentiable variable. if that is Tin, then qg is not unknown and in your simulink model it is evident that it is dependent on Tin. if the other variables are dependant on Tin, it doesnt matter. just use fcn to define the dependancy on Tin. your model is correct what's the problem?

[Muhammad Saqlain Haider]

Thank you for your responce sir. Much appreciated

qg is output of mode and actually qg is not know but I make it known by thinking that if operating set point of Tin is given then I subtracted this set point from Tin value that is coming after integral then I multiplied it with Rho_q,Va and Ca and divided it with time to make its unit in watts (as is evident from simulink model attached). Now I do not know whether its a rigt way to do or not?

[hosein Javan]

hello again. you're welcome.

your method for ode simulation is correct. in general when you want to solve "xdot=f(x,u)", you must first calculate f(x,u) then using an integrator, feed it back to your f(x,u). unless your formulation is miscalculated, the rest does not have a problem.

[Muhammad Saqlain Haider]

Thank you sir for your reply. Actually I did the same as you are saying but unfortunatelly it is creating problem in simulation. My simulations runs for some time after that it terminate the simulation and says that defirentail is finite at that time...I do not know what to do now.? any idea about that problem?.It will be much appreciated

[hosein Javan]

don't mention it. this part is not relevant to your equation. is this another equation? why is time constant "time=1"? I can't say anything for this section. you need to check the other part of the simulation where other quantities like q_lamp, q_cond ,... are defined.

[Muhammad Saqlain Haider]

This part is to calculate qg, otherwise I do not have any other option to calculate qg. Time constant is just to balance units (which should be watts for qg). Anotherthing that I need to ask, whether should I use continous time integral or discrete time integral? anyways it is giving me error in both cases but what is the rational approach according to this equation.

[hosein Javan]

you are implementing this equation to calculate qg.

I assume you want to control Tin by varying qg. and the proportional factor is:

be sure that it is correct.

about discrete integrator, no it does not make any difference normally.

[Muhammad Saqlain Haider]

Thats correct I am calculating qg like this way to control Tin. You got exactly right. I am getting this error (attached). Do you have any suggestions to resolve this?

[Sam Chak]

The model can be written in a compact form as:

If is designed as follows

then the closed-loop system becomes

Thus, at equilibrium, ,

[hosein Javan]

can you show an image of outside of this subsystem? the problem might be because of that.

the system is not stable. use a PI controller instead of rho_a*V_a*C_a after temperature error. double click on PID controller and hit "Tune". it will provide you kP and kI automatically,

[hosein Javan]

Sam Chak . The math is correct but in control case, Q_disturbance might not be directly measurable or estimated. it is just a model I suppose.

[Sam Chak]

Thanks hosein Javan. I understand that. I boldly proposed that to Muhammad Saqlain Haider because he claimed that the six q terms, which are lumped as , are known (which I interpreted it as measurable). The ideal case. Otherwise, in the non-ideal case, if is bounded and not measurable, then the PI controller as you suggested should do the job.

[hosein Javan]

Sam Chak thanks Mr.Chak. I see.

[Muhammad Saqlain Haider]

Thank you @Sam chak and @hosein Javan brothers for your support and reply. Actually hosein is right all the terms that you are calling Qdisturbance are bounded with Tin thats is why they can not be measured exclusively.

[Muhammad Saqlain Haider]

hosein Javan brother, I am not sure about PID controll. Can you eleborate it little bit more like how does it will work? As you mentioned that I should be replaced with rho_a*V_a*C_a.?

[Muhammad Saqlain Haider]

These are some errors that are existing. I guess all of them are dependent on Tin, so the real problem is in Tin. That is my understanding but the problem can be in some other thing also.

[hosein Javan]

don't mention it Mr.Haider. unfortunately I don't know from these picures. can you upload your .slx file?

[Muhammad Saqlain Haider]

Here is .slx file

[hosein Javan]

of course, I will check it tommorow, best wishes.

[Muhammad Saqlain Haider]

Ok Thank you brother for your help

[Muhammad Saqlain Haider]

Hosein Javan any update please?

[hosein Javan]

sorry I was busy today. I'll let you know a few hours later.

[Muhammad Saqlain Haider]

Sure brother..Appreciate it

[hosein Javan]

sorry, I'm using matlab 2016a. can you export it to this older version?

[Muhammad Saqlain Haider]

@hosein Javan I hope it will work now

[hosein Javan]

yes it worked now. the model is too complicated to comprehend, I suggest first try to simplify the model in this way:

the model can be represented by "n" ode in the form of "Xdot = f(X,U)". so create a one subsystem that has input "X,U" and outputs "f(X,U)". then feed it to an integrator and then feed it back to "X".

also, some "from workspace" blocks are undefined. I suggest instead of these, use "constant" with pre-defined values. I saw that some of them were tout. how is it that tout is needed for a simulation? use "clock" if necessary however in most cases it is unneeded.

Once you know that everything that is related to model is done, create a subsystem that defines your system. this subsystem has input "U" and outputs "X".

right now the model is ready for simulation. provide some arbitrary inputs for "U" and check if the system response "X" is expected, and no error like infinite states apper. note that the "X" must be stable at equilibrium. if not the system is ill-defined. for now just the model is defined. please don't add any extra feedback closed-loop controllers, and only check its response.

allrigth, when system check done, it is time for control, but for now, check only the system. I suggest you once again write all equations on paper and create a new file based on above steps. once you are done with all that please let me know.

[Muhammad Saqlain Haider]

As for as this is concerned "the model can be represented by "n" ode in the form of "Xdot = f(X,U)". so create a one subsystem that has input "X,U" and outputs "f(X,U)". then feed it to an integrator and then feed it back to "X".".

I am doing same, Tin is my output of model and I am feeding it back because in all inputs it is needed. As for as ode is concerned I do not think so there is any error in writing equation in simulink but yes there is something wrong because of which model is not working

[Muhammad Saqlain Haider]

[hosein Javan]

yes I know. I just meant that the model is confusing for a lot of line signals. that's why I suggested at first try to simplify the model. you can also use "Goto" and "From" blocks to reduce its complexity. please test the model once without the feedback controller to ensure that the error is not due to system definition. good luck.

[Muhammad Saqlain Haider]

Thank you so much brorher for your support. I found one problem that is very obvious, it is because of division by zero error. Do you have any idea how to fix this problem because itis happenening in some equations? @hoseinJavan

[hosein Javan]

don't mention it my brother. division by zeros normally should not happen. probably some algebraic equations are defined in such a way that some variable is divided by a state variable at initial condition which is zero. please first check that this signal is defined correctly, if it is correct try to:

1. if there is way to define this equation in another way that does not require division by zero, do it. for example if the equation is y1=x2/x1, and we already know that x1=1/x3, then we can rewrite y1=x2*x3.
2. some algebric equations have a limit at zero. for example limit(sin(x)/x)=1. if there is any equations that have limit, try to add a switch that when the signal is zero, produce this pre-calculated limit.
3. change initial condition in a way that division by zeros does not happen. for example if y1(t)=x2(t)/x1(t) and x1(0)~=0, then division by zeros does not happen.
4. if none of these methods are incapable, we can use saturate the value of denominator dynamically to "eps", this way a very large but not inf occurs. however it is the least we can do and it is strongly recommended to avoid this method, because it might lead to unwanted phenomena like many zero-crossings and instability of the system due to very large signal.

[Muhammad Saqlain Haider]

hosein javan ,Great thank you so much brother your answers are helping me alot. As we know that Simulink runs simulation in seconds but this is a model which I need to run for hours and days. I wanted to give parameters values on hourly basis in my model as is given in my data but I know that Simulink runs simulation in seconds so any suggestions that How can I give these values in simulink input. I am attaching the image to show you the nature of data. Time here is in hours hence every value

[hosein Javan]

suppose that simuliink measures time in terms of seconds. then it means that if state variable x(t) is in terms of meter then the state variable v(t)=dx(t)/dt is in terms of meter per seconds at any seconds.

now imagine that we don't even know that simulink measures time in seconds, for example we might think that simulink measures in hours. then if state x(t) is distance in terms of meter, this time, v(t)=x(t)/t is speed in terms of meter/hours.

so we see that only units that are time-based change by a scale. now for example here wind speed is in terms of a distance unit like meter or kilometer or miles,... per hours.

so in a conclusion, if we suppose time is in terms of hours, only the time depandant units change by a factor that can be easily retreived, the other units don't change.

[Muhammad Saqlain Haider]

So, it means that if I am running simulation for 24 units. Should I consider it as 24 hours?. Will it fullfill the purpose?

[Muhammad Saqlain Haider]

One more this, for example I put stop time as 24 and considering that its in hours so it means that model is running for 24 hours means for one day.Now if I wanna run my model for couple of days what shoud I do? Like for 2 days should I put 48 as a stop time? By considering that I am giving time in hours? @hosein Javan

[hosein Javan]

Yes, Exactly. Just turn every thing that you thought was seconds to hours. Only units that deal with time like speed. But the other units that are not a "time-rate of something" like heat and temperature are no time-depandant units. Remember that some quantities have time hidden in them. like electric current which is time rate of charge flowing.

[Muhammad Saqlain Haider]

Thank you so much brother @hosein Javan. I really appreciate your quick and wonderful responce.

[hosein Javan]

don't mention it brother. please if you are convinced with the answer, hit "accept".

[Muhammad Saqlain Haider]

How can I accept the answers? I did not find any tab here. Can you please guide me to do so? @hoein Javan. Another question is related to PID controller, can you please guide me after looking at model that where can I setup a PID control and tune it in my model? Can you please eleborate me through picture that how can I setup? I saw the videos but I did not get idea because they did not use set point but in my case I have set point. It would be really appreciated as well. Thank you

[Sam Chak]

Hi Muhammad Saqlain Haider,

Have you successfully run and stabilize the system with the designed controller without any singularity (division-by-zero) error message?

[Muhammad Saqlain Haider]

@Sam Chak, not yet. Actually I am trying to built a PID controller now but I am not sure how can I setup in my model. I do not know what will be my disturbance because I know setpoint which is 20 but I also need disturbance for PID. I am sending a picture here from where you can see why I need disturbance. Guide me if you can please

[Sam Chak]

@Muhammad Saqlain Haider. The disturbance blocks in the images are "External disturbances", are external elements in the process that perturb the responses of the system states. External disturbance can be a constant or time-varying, or a mixed of both in different periods of time.

External disturbance such as gravity (gravitational force) is generally a constant (for land dwellers) and is independent of the system states.

External disturbance such as air drag (aerodynamic force due interaction with external air particles) depends on the system's velocity (one of the states in a motion system).

[Muhammad Saqlain Haider]

Thanks for ypour responce Sir. What could be the disturbance in my case? I am attaching the image of my model so Can you give me aqny idea that what are some disturbances in my model and secondly what does transfer function means in my model? Does it means integral that I am using including all the terms before integral? Thats what I understand but i can be wrong

[Sam Chak]

Hi Muhammad Saqlain Haider, I didn't study your system dynamics in the Simulink model. Perhaps Dr. hosein Javan can provide constructive insights on this matter because he has previously viewed and commented on your Simulink slx file.

Answer 1

[hosein Javan]

hello again. this should look like this:

by clicking on "Tune" an app opens and tries to first linearize the system and then choose the best coefficients.

[Muhammad Saqlain Haider]

Great Thank you Sir @hosein Javan. Detailed reply as always. May I ask that what are meaning of linearization of system?

[Sam Chak]

Hi Muhammad Saqlain Haider, please note that auto-tuning works if there are no discontinuity or singularity issue in the system process. Of course, the tuning results are only meaningful if you correctly place the PID controller block in the control loop.

[Muhammad Saqlain Haider]

Yeah right sir, that is what I do not know where I have to put PID controller?. As I am I have to control Tin (inside temperature) by changing qg so I do not know how and where I have to put PID. As Dr. @hoesin javan mentioned in above picture, it is not working actually and I think it should not work because it is not controlling inside temperature which is Tin and there is no feedback loop as well. So I am not sure, if you have any idea please help

[hosein Javan]

you're welcome Mr. Haider. many systems are nonlinear, that means they are not represantable by a linear differential equation. PID controller is only for linear systems, so in order to use that for a nonlinear system, the app tries to linearize the system first. meaning that for a small perturbation around the equilibrium, it is approximated by a linear system. as Mr. Chak said, the block tuning will work in certain condition, I might add that there could be other reasons why this fails like high nonlinearity and ofcourse the system ill definition. the PID is placed in the right place where after T compares with setpoint value and produces the qg necessary to acquire the wanted temperature. if tuning did not work, try to retune it yourself. first change the kP factor with having kI=0. remeber that kP certainly does not make your quantity to reach its setpoint but it should gets near. after that try to gradually increase kI, but before all these, I suggest that run the simulation once without controller, and try to manually change qg and see that the temperature reqires how much qg to estimate how much the controller should produce.