Approximating Differential Equation on the Basis of Measured Values

I'm back. OK JCH, one more chance. This should be simple. The control output is 1 and there are only 26 points just how you like it. The data is in CSV, comma separated variable, file format.

Time, Position, Control Output

0,0,1 0.1,0.0280961296169728,1 0.2,0.105480069053459,1 0.3,0.22321745414104,1 0.4,0.373993446175833,1 0.5,0.551819161757164,1 0.6,0.751791317868303,1 0.7,0.96989544591241,1 0.8,1.20284477699198,1 0.9,1.44794833233238,1 1,1.70300292485492,1 1.1,1.9662047375435,1 1.2,2.23607692993412,1 1.3,2.5114103673215,1 1.4,2.79121509393783,1 1.5,3.0746806025518,1 1.6,3.36114330596755,1 1.7,3.65005990494049,1 1.8,3.94098558367094,1 1.9,4.23355615778425,1 2,4.5274734583331,1 2.1,4.82249336523072,1 2.2,5.1184160098546,1 2.3,5.41507775361695,1 2.4,5.71234462057353,1 2.5,6.01010692049863,1

Have fun.

Peter Nachtwey

Reply to
Peter Nachtwey
Loading thread data ...

"Peter Nachtwey" schrieb im Newsbeitrag news:S4mdnYNVEokMvqDanZ2dnUVZ snipped-for-privacy@comcast.com...

See Page 1/2

formatting link

Reply to
JCH

Newsbeitragnews:S4mdnYNVEokMvqDanZ2dnUVZ snipped-for-privacy@comcast.com...

It isn't right. Not even close. Did you plot the estimated PV against the actual PV? Then it would be obvious.

I will post the solution after the week end if you don't find the solution before then.

Peter Nachtwey

Reply to
pnachtwey

I did it again without normalization. Basically the same solution:

Time span is now real 0 to 100:

formatting link

Reply to
JCH

Newsbeitragnews: snipped-for-privacy@s6g2000prc.googlegroups.com...

100:
formatting link
It is still wrong. Actually, I can't follow what you are doing. I don't know why you show disturbances because there are none. I don't know why you have PID gains when the model is wrong. My data is an open loop response to a step control output of 1. There is no noise, no dead time, no offset and only one time constant (hint). It is that simple. You didn't plot your estimated response and the given response to compare them.

Peter Nachtwey

Reply to
pnachtwey

Newsbeitragnews: snipped-for-privacy@s6g2000prc.googlegroups.com...

If no disturbance then there is nothing to control. I show a close loop control.

That is just PI-controlled with only 1 disturbance (z2):

formatting link
This is feedforward PID-controlled with 2 disturbances (z1, z2):
formatting link
The estimated plot of process transfer function is black colored. The response is red colored but is covered by the black line.

Reply to
JCH

- snipped-for-privacy@s12g2000prg.googlegroups.com...

roups.com...

Your model, 0,1796511 y ' ' + 0,3153461 y ' + 0,01438058 y =3D 1, is still wrong. It is not even close. At steady state y' and y'' will be

0 there for y =3D 1/0,01438058 so y =3D 69.5. What is y? Is it a position or a velocity? If y is a position I can tell you it is wrong because as long as the control signal of 1 volt is applied the positoin will keep changing at not stop at 69.5. There is no steady state value for positon. If y is a velocity then you are indicating that the steady state velocity is 69.5 position units/second. This is clearly not right because you an the the velocity approaches steady state value of 3 poisition units per second from my data. Look at the last two data points. 2.4,5.71234462057353,1 2.5,6.01010692049863,1

The difference in time between the two points is 0.1 seconds and the difference between the two points is 6.01010692049863-5.71234462057353 is about 0.298 so the velocity is 2.98 positoin units per second.

I have now given away what the answer is. Truthfully, I could get a good estimate of what the transfer function is by just looking at a graph of the data this model is that simple.

You still haven't graphed the estimated position versus the actual position data given to see if your model matches the actual position data. If the model isn't right then the tuning values aren't right. Concentrate on the system model.

Peter Nachtwey

Reply to
pnachtwey

Newsbeitragnews: snipped-for-privacy@s12g2000prg.googlegroups.com...

Newsbeitragnews: snipped-for-privacy@s6g2000prc.googlegroups.com...

Your model, 0,1796511 y ' ' + 0,3153461 y ' + 0,01438058 y = 1, is still wrong. It is not even close. At steady state y' and y'' will be

0 there for y = 1/0,01438058 so y = 69.5. What is y? Is it a position or a velocity? If y is a position I can tell you it is wrong because as long as the control signal of 1 volt is applied the positoin will keep changing at not stop at 69.5. There is no steady state value for positon. If y is a velocity then you are indicating that the steady state velocity is 69.5 position units/second. This is clearly not right because you an the the velocity approaches steady state value of 3 poisition units per second from my data. Look at the last two data points. 2.4,5.71234462057353,1 2.5,6.01010692049863,1

The difference in time between the two points is 0.1 seconds and the difference between the two points is 6.01010692049863-5.71234462057353 is about 0.298 so the velocity is 2.98 positoin units per second.

I have now given away what the answer is. Truthfully, I could get a good estimate of what the transfer function is by just looking at a graph of the data this model is that simple.

You still haven't graphed the estimated position versus the actual position data given to see if your model matches the actual position data. If the model isn't right then the tuning values aren't right. Concentrate on the system model.

See the 2.2*10^-4 % Model: New

formatting link
Compare with the 1.4 % Model: Old
formatting link
Can you see any difference in control performance?

Reply to
JCH

- snipped-for-privacy@e10g2000prf.googlegroups.com...

roups.com...

egroups.com...

New 0,1666572 y ' ' + 0,3333362 y ' + 2,174717E-06 y =3D 1

this is still wrong.

Why do you worry about controlling the wrong model? The controller and the model are wrong. Stick to the problem. Find the system gain and time constant. Let common sense guide the math and not the other way around.

Another hint. Differentiate my data. Then you will have velocity as a function of time. I am pretty sure you can find the solution then.

Peter Nachtwey

Reply to
pnachtwey

Newsbeitragnews: snipped-for-privacy@e10g2000prf.googlegroups.com...

Newsbeitragnews: snipped-for-privacy@s12g2000prg.googlegroups.com...

Newsbeitragnews: snipped-for-privacy@s6g2000prc.googlegroups.com...

New 0,1666572 y ' ' + 0,3333362 y ' + 2,174717E-06 y = 1

this is still wrong.

JCH

I regard a steady state error of ~2*10^-4 % as more than accurate for technical applications.

SHOW YOUR ANSWER!

Peter

Another hint. Differentiate my data. Then you will have velocity as a function of time. I am pretty sure you can find the solution then.

JCH

I see the differentiation whenever I look at my program.

See velocity: Page 2

formatting link

Reply to
JCH

Newsbeitragnews:858a500f-3f26-437e-8506-

Error between what and what? What is steady state? I have asked before, what is y? Is it a position or a velocity. In either case your equation is wrong. From you formula

0,1666572 y ' ' + 0,3333362 y ' + 2,174717E-06 y = 1

the steady state value of y, what ever that is, is 1/2,174717E-06 which is a very big number. It is way too big for a steady state velocity and there is no steady state position. I don't see how you expect that to match the given data. I bet you still haven't graphed your model and compared it with the given data. Do you want me to do it for you?

Peter Nachtwey

Reply to
pnachtwey

The answer is, approximately:

H(s) =3D 8.94 ---------------------------- s^2 + 4.904 s + 8.94

That was fun.

Reply to
Andrew

I am not sure if I understand what is being asked here. Are we trying to find the equation which produced that set of data? Or are we trying to find a differential equation for which that set of data is the solution?

I'll assume that we are looking for the equation which defines this set of data:

In the s domain, this set of data can be approximated by the following:

H(s) =3D 8.94 ------------------------- s^2 + 4.91 + 8.94

In the time domain, this is equivalent to:

H(t) =3D 1 - 1.75*e^(2.45t)*cos(1.71t-55)

Is this what is being asked for here?

Andy

Reply to
Andrew

Good question. The goal is to find the transfer function that best estimates plant's response. Normally one captures the process variable, PV, as a function of the control output, CO, and time. In JCH's data the CO is a step jump of one.

If you use AB PLC's then you can save the PV as a function of CO in a trend file.

So why bother? If you know the plant transfer function then one can calculate the PID gains that will provide almost any desired response. Calculating the PID gains is a different issue but the difficult part is finding the plant transfer function. After that calculating the PID gains is easy. See

formatting link
for the formulas for calculating the gains.

Andrew, I don't see how you calculated H(t). It doesn't look right, but the transfer function

is pretty good for a two pole solution. See this ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20SysID%20JCH%20Andrew.pdf You can see that your solution is different from mine but that could be because you used a different evaluation function. I basically use the sum of squared errors or divide that by the number of points to get a mean sum of squared error..

You should know that JCH posted actual transfer funciton back on Oct

  1. The actual plant is 5 repeated poles with a time constant of
0.123.

Are you familiar with Runge Kutta and optimizing functions? So far we are only at the beginning stages. Issues that haven't been covered yet include offsets in the plant transfer function, using foricing functions other than 1, and dead time.

Peter Nachtwey

Reply to
pnachtwey

I am not sure if I understand what is being asked here. Are we trying to find the equation which produced that set of data? Or are we trying to find a differential equation for which that set of data is the solution?

I'll assume that we are looking for the equation which defines this set of data:

In the s domain, this set of data can be approximated by the following:

H(s) = 8.94 ------------------------- s^2 + 4.91 + 8.94

In the time domain, this is equivalent to:

H(t) = 1 - 1.75*e^(2.45t)*cos(1.71t-55)

Is this what is being asked for here?

Your's equivalent to

1 H(s) = ----------------------- 0.112*s^2 + 0.549*s + 1

My approximation for 2nd order

1 H(s) = ----------------------- 0.103*s^2 + 0.421*s + 1

That is similar to your solution.

See difference to 5th order:

formatting link

Reply to
JCH

snipped-for-privacy@e10g2000prf.googlegroups.com...

I looked up the 5th order equation. It is clearly a perfect fit. My second order solution was easy. All you need to define a second order system is the natural frequency and the damping ratio. You can find that with a small amount of analysis from the data. But I also assumed that it was a 2nd order system. I have a couple of questions?

  1. How would you know it was a 5th order system?
  2. How do you find the parameters for a 5th order system?
  3. What is the time domain response?
  4. Other than a good math problem (which is always a good thing), why would you model a system as 5th order when a 2nd order approximation is apparently close enough?

andy

Reply to
Andrew

om/public/NG/Mathcad%20-%20SysID%20JCH%20Andr...

Peter, The relationship between H(s) and H(t) is a direct translation from the laplace domain to the time domain. If you plot them together, you will see that they are identical.

Andy

Reply to
Andrew

snipped-for-privacy@e10g2000prf.googlegroups.com...

I didn't at first. JCH said that we should limit the number of poles to 3. First I plotted JCH's data and I can tell by the inflection that is was at least 3 or more probably 5. It takes only a few seconds to try many options and compare the error. The model with the smallest error between the actual and estimated data is usually the right one.

All the poles were at the same location so it was easy. I use a optimizer function that minimizes the error between the actual and estimated data. In the link below this function is called Minerr. Scilab also has lsqrsolve and optim that do the same things.

I posted a link to my solution on Oct 22. ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20SysID%20JCH.pdf notice it took little effort to change n and find a new solution. When the error stops getting smaller what is left is noise and un- modeled non-linear features. The equation for the time domain response follows a pattern. I don't normally use equations in the time domain though. I usually use a state space method or RK4 to integrate the differential equation(s). Note the (s).

That is a good question too. Is a second order approximation good enough? I think this question deserves its own thread though. A second order approximation may be good in some cases but not others. This is an opinion, not fact and the answer may depend on the definition of what is good enough. I don't think, but haven't shown, a second order approximation is good enough in JCH's example because all the poles are at the same place. There are no dominate poles with the other three being high frequency poles that will barely affect the system. In any case wouldn't you want to know were the poles are before you decide to ignore them? I know that I can design a better controller for JCH's 5 pole system knowing it is a 5 pole system than if I assume there are only two or three poles. Even if you only have a PLC with a PID I still think you are better off knowing reality

I use system identification in motion control. Most system can be modeled as a type 1 first or second order system. Try my test data. Try the easy on first with only 26 points. Both sets are more typical of what one sees in a motion control application than JCH's test data.

It is about time some the lurker come out. I can tell by the down loads from my ftp site that this thread isn't being ignored.

Peter Nachtwey

Reply to
pnachtwey

Newsbeitragnews: snipped-for-privacy@e10g2000prf.googlegroups.com...

I looked up the 5th order equation. It is clearly a perfect fit. My second order solution was easy. All you need to define a second order system is the natural frequency and the damping ratio. You can find that with a small amount of analysis from the data. But I also assumed that it was a 2nd order system. I have a couple of questions?

  1. How would you know it was a 5th order system?
  2. How do you find the parameters for a 5th order system?
  3. What is the time domain response?
  4. Other than a good math problem (which is always a good thing), why would you model a system as 5th order when a 2nd order approximation is apparently close enough?
Reply to
JCH

Newsbeitragnews: snipped-for-privacy@i12g2000prf.googlegroups.com...

ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20Type%201%20First%20Order%20System.pdf

One can see from the graph that the velocity increases and is reaching a steady 3 position units per second. Therefore one can see the gain is about 3 positions per second per control output unit. If one plots the velocity as a function of time one can see this is a first order system because there is no inflection point. in the response. One can also see that it takes about 2.5 seconds to get within 1% of 3 position units per second so the time constant must be around 0.5 seconds. Before I could do all the math I relied heavily on just analysing the graphical data.

One can see there is a big limitation to being limited to analysing only 25-30 points. One can see that the response to a step change takes about 5 time constants so that leaves only 5 to 6 samples per time constant. That isn't a high enough sample rate.

There are still some issues not addressed. This data is perfect. What if it one added some noise and quantized the data to reflect converting the feedback in to a digital format. What happens if there is a offset? In open loop system will often drift even when the output is about 0. This means one must add another term so the differential equation becomes

vel'=-alpha*vel+Gain*alpha*u(t)+offset

One can not identify a system like this when u(t) is just a step of

  1. u(t) needs to change between different levels so the system identification can tell the difference between the output due to u(t) and the offset. No system identification can be done on a real system with a offset if the system is only excited with one step change.

Peter Nachtwey

Reply to
pnachtwey

PolyTech Forum website is not affiliated with any of the manufacturers or service providers discussed here. All logos and trade names are the property of their respective owners.