Approximating Differential Equation on the Basis of Measured Values

Approximating Differential Equation on the Basis of Measured Values, not physical derivation.
When starting up technical systems one can utilize a plant step transfer
function for getting data like (x = time, decimal commas)
x y
0 0 0,087 0,001 0,174 0,011 0,262 0,066 0,349 0,158 0,436 0,289 0,524 0,415 0,611 0,560 0,698 0,661 0,786 0,764 0,873 0,834 0,961 0,895 1,048 0,922 1,135 0,956 1,223 0,968 1,310 0,988 1,397 0,978 1,485 0,998 1,572 0,994 1,660 1,004 1,747 0,994 1,834 1,003 1,922 0,994 2,009 1,001 2,096 0,996 2,184 1,000 2,271 0,994 2,359 1,001 2,446 0,990 2,533 1,000
I searched the internet for getting the appropriate differential equation, e.g.
F(s)=(T*s +1)^n
or
A3*y''' + A2*y'' + A1*y' + y = 1
and could'nt find a website.
T=? n=?
or
A1=? A2=? A3=?
Can someone provide the solution for comparison under the condition
n<=3
?
TIA
--
Regards/Gre http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann eMail aktuell: snipped-for-privacy@nospam.arcornews.de
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
This looks more like a test since we know you have a program to solve these kinds of problems. Good test though. I will work on it when I have time. It will be interesting if anyone else will be up to the challenge.
Peter Nachtwey
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
On Mon, 22 Oct 2007 13:45:44 +0200, JCH wrote:

You've searched on the key words "system identification"?
--
Tim Wescott
Control systems and communications consulting
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

Can't you see? It is a test for all to try to solve. If you have looked at JCH's website you can see that he has a program that will do system identification. Yes, we are going to be playing JCH's game on his field. However, it is a worthy problem because if you can master system identification then tuning takes much less tweaking and coffee drinking. Game on.
Peter Nachtwey
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
On Mon, 22 Oct 2007 09:27:47 -0700, pnachtwey wrote:

I don't wish to get drawn into a 'game', so when he's being reasonable I answer his posts for the edification of others if not him, and when he gets all weird I ignore them. So far he hasn't gotten bad enough for me to plonk him, but that may come.
Game _off_.
--
Tim Wescott
Control systems and communications consulting
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

JCH is being reasonable. He has provided enough information to get an answer and his test data is pretty simple. Here is my partial solution: ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20SysID%20JCH.pdf I think it is important that people give this a try. System identification is one of those things that separates the tweakers and coffee drinkers from the real control engineers. I think this is important enough that I will show how I calculate my solution sometime after the next weekend. I don't want to ruin the experience of the aha discovery experience for anyone.
Note my solution assumes that all the poles are specified by just one time constant. One may be able to get a better solution if one tries a different time constant for each pole.
Peter Nachtwey
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
pnachtwey wrote:

So riddle me _this_, Peter.
If you are tuning a process that's nonlinear, and whose apparent linear model therefore necessarily changes with tuning, wouldn't it be necessary to repeat ones tests as one changes ones tuning?
And if one were smart enough to repeat ones tests to get a refined model, and refine ones tunings as a consequence of the apparent change in the plant model, wouldn't it be intellectually honest to recommend just that procedure to others?
And if this process takes time, wouldn't it be wise and kind to recommend patience in the taking of the data?
And if, in advocating patience in taking one's data, one mentioned that one needed to tweak (or in other words to re-tune based on one's refined model) and drink coffee (or in other words to not rush to an incorrect tuning based on a too-hurried measurement process), wouldn't one expect true professionals to treat this advice with respect, even if they somehow disagreed that one should engineer from complete and correct data?
And then, would one consider endless derogatory mentions of "tweaking and coffee drinking" from one of ones peers to be a bit, well, _un_professional?
Please explain _that_, while you're playing games with JCH.
Thanks in advance.
--

Tim Wescott
Wescott Design Services
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
Tim Wescott wrote:

And thank you, Tim, For so eloquently saying what I was thinking. Mathematical solutions are starting points that are suspect till proven effective by actual trials.
A control engineer who often tweaks and drinks tea, refining an approximate solutions (approximate because of linear assumptions) on nonlinear systems,
John Popelish
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
wrote:

See comparison http://home.arcor.de/janch/janch/_news/20071023-comparison /
Peter, your result: Fig.1
F(s)=(0.206977*s+1)^3
equivalent to
0.01042659*y''' + 0.1323628*y'' + 0.6381982*y' + y = 1.00
This is a valid but not a good result.
A better result: Fig.2
--
Regards/Gre http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann eMail aktuell: snipped-for-privacy@nospam.arcornews.de
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

I question your numbers

My Mathcad computes. I just cut and pasted these number from my Mathcad.
.886679e-2*s^3+.128518*s^2+.620931*s+1. This is more than simple round off error. I agree the results for n=3 aren't very good. I think it is strange that your result for n=3 looks like my n=5. You should recheck how you calculate your numbers. 0.206977^3=0.008867 not 0.0104...
Peter Nachtwey
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
wrote:

That's your 3rd oder result. See above.

Solutions:
Peter: 0.008867000*y''' + 0.128518*y'' + 0.620931*y' + y = 1.000 JCH : 0.009598422*y''' + 0.104815*y'' + 0.541117*y' + y = 1.000
That's the first time I show my numbers.

In addition the program finds boundary values (start values y''(0), y'(0), y(0)) so that the result in Fig.2 can be found as is shown: http://home.arcor.de/janch/janch/_news/20071023-comparison /
These boundary values should be used in the target filter F3(s) preventing overshooting u: http://home.arcor.de/janch/janch/_control/20070613-pd2 (pid)z1z2/
The basic idea is best control for PT5 system.
--
Regards/Gre http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann eMail aktuell: snipped-for-privacy@nospam.arcornews.de
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

Your solution has two imaginary poles. It is not in the form of 1/(T*s +1)^n as indicated on your website. My solution current solution is limited to real poles as on your website. I see what I come up with when I allow imaginary poles too.
Peter Nachtwey
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
wrote:

Ok, therefore I used ???, not T.

0.009578639 y ' ' ' + 0.1046687 y ' ' + 0.5406281 y ' + y = 1.000227
Start values
y (0) = -0.000878999 y ' (0) = 0.2983399 y ' ' (0) = -10.97643
Solve this ODE and you have the solution/approximation shown Page 1 (Fig.2) and Page 2 in:
http://home.arcor.de/janch/janch/_news/20071024-comparison /
--
Regards/Gre http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann eMail aktuell: snipped-for-privacy@nospam.arcornews.de
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

STOP RIGHT THERE!!!! The forcing function and the initial conditions are not part of the model. You can't optimize them just to make your model look better for just this set of data. The model should work for any set of data and initial conditions. ( We are assuming the plant is linear now ). The goal is to find the best model, not make the model response follow the actual response at the cost of everything else. In this case it is reality. You must start your system at a quiescent state or know what the initial state is but the initial state is NOT a part of the parameters that is optimized. The same goes for the forcing function. Your initial problem set the forcing function to 1 not 1.000227. The forcing function is an input and not part of the model and should be left at one. The initial conditions should be left at 0. You were doing well. Don't go weird on us now.
Peter Nachtwey
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
wrote:

In a 'context of getting better results' one must use them. It's up to you to ignore them.
They are numbers just different from zeros. I can easily simulate the usage on my computer. That means move the target value (u) using initial values, and you move it like having an ODE of order 5.
Example: built-in a target filter F3(s) 3rd order but acting like 5th order: http://home.arcor.de/janch/janch/_control/20070613-pd2 (pid)z1z2/
--
Regards/Gre http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann eMail aktuell: snipped-for-privacy@nospam.arcornews.de
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

One does not optimize the inital conditions. When you do you are playing mathematical gains but you aren't doing system identification. I noticed you made no comment about optimizing the forcing function. Optimizing the the inital conditions makes as much sense as optimizing the forcing function.

To you this is just a numbers game but in reality it does no good if you aren't learning the about the system. Personally, I would use the 5 repeat poles model. I don't see how using a wrong model will be beneficial do designing a controller.
Peter Nachtwey
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
[...]

The step was chosen 1. That's all. The approximation routines changed it a little bit. Why bother about 0.000227=0.00000227% difference?

Look again at the Fig.2, page 1 and Fig. at page 2: http://home.arcor.de/janch/janch/_news/20071024-comparison /
My alternate '5 repeat poles model' approximation: 1/(0.122966+1)^5 2.811421E-05y'''''+0.00114317y''''+0.01859326y'''+0.1512065y''+0.6148302y'+1=1
That's what you can easily find but is more difficult to handle in reality, I guess.
--
Regards/Gre http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann eMail aktuell: snipped-for-privacy@nospam.arcornews.de
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
Newsbeitrag

That is good because our solutions are very close if you really meant: 1/(0.122966*s+1)^5 (I added the s ) and 2.811421E-05y'''''+0.00114317y''''+0.01859326y'''+0.1512065y''+0.6148302y'+ y=1 ( I change the 1 to a y )
Hopefully you didn't try to optimize the initial conditions. So: 1 what equation did you use to generate the test data? 2. What were the initial conditions for you generated test data? 3 Doesn't your optimize algorithm return a parameter the indicates how good the fit is like a norm, sum of squared errors or mean squared error? This term can be handy because it is the high bound of the covariance needed for setting up the Q matrix of Kalman or H infinity filters. The system identification MSE is the combination of the process and measurement errors.
I am using Mathcad. It has a built in optimizer called minerr that I use to quickly do system identification and other optimization problems. I believe it use the BFGS algorithm. In Scilab I use lsqrsolve which uses the Levenberg-Marquardt algorithm. Scilab also has a optimize function Which optimizing algorithm do you use?

Yes, but how hard is it for you do generate feed forward gains for the 5th derivative?. From what I have seen it wouldn't be hard. How hard is it to generate a motion profile where the 5th derivative is smooth at the beginning and ending endpoints for a ramp? I don't think you would have much problem with that either. Generating continuous and smooth 5th derivatives can be done but it would definitely be a custom job. I think it requires a 9th order polynomial. In the plot on your website the PID2D wasn't necessary most of the time as you relied almost entirely on the target filter and feed forwards to generate the control output. You would do a similar thing here except you should drop your ramp filter in favor of a high order polynomial that will allow the 5th derivative to be 0 at the end points.
This means that: y(0)=StartPoint, y'(0)=0, y''(0)=0, y'''(0)=0, y''''(0)=0, y'''''(0)=0 y(t)=EndPoint, y'(t)=0, y''(t)=0, y'''(t)=0, y''''(t)=0, y'''''(t)=0 t = ramp time
I am sure you can do that. The point is that a smooth ramp and a simple PID may be enough because almost all the control output will be from the feed forward values. The PID only needs to correct for modeling errors which are small and disturbances which there are none in your test case. Not all the poles can be placed but a few can be. These can be made dominate. The rest will just float but if they aren't excited then they should cause little harm. First calculate the ramp and the feed forwards.
When generating ramps or motion profiles, I have never had the need to do more than keep the 4th derivatives continuous and the 3rd derivative smooth in actual practice. Most people think this is extreme. Most motion controllers use only 3rd order motion segments to piece together a motion profile. 3rd order segments are only smooth at the derivative or velocity and continuous at the second derivative or acceleration.
Peter Nachtwey
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

Sorry for the editing errors.

1/(0.123*s+1)^5 (+/-1% error for making real)

Initial values zero.

Sum(|errors|) => MIN, automatically, but still manually adjustible if necessary

... easily on a computer or computer-like equipment (industrial PCs).
In my mind I wanted to limit it to standardized 3rd order equipment without loosing control quality. Therefore one have to calculate the initial values for using an 3rd order ODE equivalent. This is the point I was interested in.
See also other people looking for start values: http://www3.imperial.ac.uk/portal/pls/portallive/docs/1/7290459.PDF
- The Reduced-Order Observer Using the Application of Matrix Generalized Inverses: Page 15
- Finding the Optimal Initial Condition for Time-Varying Systems by Using the Min Max Optimization: Page 19
In earlier discussions I didn't mention observers. But mathematical models are observers. I used a State Observer. Necessary correcting is done by a PID-controller.
Additional comparing real PV with mathematical model PV and proportional correcting the mathematical model will be a Luenberger observer.
Addional information: http://home.arcor.de/janch/janch/_control/20071028-orderchange /

That's why I said in an other thread that Luenberger observer is not necessary. The advantage is neglectible

One can also utilize static and dynamic disturbance compensation: http://home.arcor.de/janch/janch/_control/20070613-pd2 (pid)z1z2/ See Page 1, z in formulae

I wouldn't recommend ramps with constant velocity. One must avoid jerks. The best results can be found using the target transfer function F3(s) identical to the process transfer function F1(s).

The most important preposition is to keep all variables in control range at any time (bumpless control). This was not accomplished in Page 6.
Just PID-control: Page 6 http://home.arcor.de/janch/janch/_control/20070613-pd2 (pid)z1z2/
--
Regards/Gre http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann eMail aktuell: snipped-for-privacy@nospam.arcornews.de
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
JCH wrote:

I think you should try squaring each error and summing it and minimizing the sum of square error ( SSE ). I like to divide by the number of points do normalize the error ( MSE ). This doesn't affect the optimizing at all. Obviously if there are more points the error will be bigger. Dividing the sum of square error by the number of points gives me a better idea of how close the fit is.
As far as making things easy. I won't be. Even your 3 poles solution is too complicated to be used on a PLC because they only have PIDs, not PID2D closed loop control. A PID can only be used to place two poles because the I gain comes with its own pole. In total there will be 6 poles but only two can be placed.
I will see if I can find a suitable control method for your 5 repeated poles. Obviously the plant is already stable so the goal would be to be able to get the control the PV from 0 to 1 faster than what can be done in open loop. I will use the my model to calculate the pole positions but run the simulation using your transfer function used to generate the test data because that is the actual system. At this time the sun is out and must be enjoyed so maybe tomorrow.
Peter Nachtwey
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

Polytechforum.com is a website by engineers for engineers. It is not affiliated with any of manufacturers or vendors discussed here. All logos and trade names are the property of their respective owners.