Approximating Differential Equation on the Basis of Measured Values

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

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
http://www.wescottdesign.com

Need to learn how to apply control theory in your embedded system?
"Applied Control Theory for Embedded Systems" by Tim Wescott
Elsevier/Newnes, http://www.wescottdesign.com/actfes/actfes.html

text -

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

On Mon, 22 Oct 2007 09:27:47 -0700, pnachtwey wrote:


text -

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
http://www.wescottdesign.com

Need to learn how to apply control theory in your embedded system?
"Applied Control Theory for Embedded Systems" by Tim Wescott
Elsevier/Newnes, http://www.wescottdesign.com/actfes/actfes.html

text -

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

pnachtwey wrote:

text -

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
http://www.wescottdesign.com

Do you need to implement control loops in software?
"Applied Control Theory for Embedded Systems" gives you just what it says.
See details at http://www.wescottdesign.com/actfes/actfes.html

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

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/Grüße    http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann  eMail aktuell: janch@nospam.arcornews.de
                 Microsoft-kompatibel/optimiert für IE7+OE7

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

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/Grüße    http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann  eMail aktuell: janch@nospam.arcornews.de
                 Microsoft-kompatibel/optimiert für IE7+OE7

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

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/Grüße    http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann  eMail aktuell: janch@nospam.arcornews.de
                 Microsoft-kompatibel/optimiert für IE7+OE7

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

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/Grüße    http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann  eMail aktuell: janch@nospam.arcornews.de
                 Microsoft-kompatibel/optimiert für IE7+OE7

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.


order:http://home.arcor.de/janch/janch/_control/20070613-pd2 (pid)z1z2/
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

[...]

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/Grüße    http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann  eMail aktuell: janch@nospam.arcornews.de
                 Microsoft-kompatibel/optimiert für IE7+OE7

2:http://home.arcor.de/janch/janch/_news/20071024-comparison/

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

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/Grüße    http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann  eMail aktuell: janch@nospam.arcornews.de
                 Microsoft-kompatibel/optimiert für IE7+OE7

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