Approximating Differential Equation on the Basis of Measured Values

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

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:

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

"pnachtwey" schrieb im Newsbeitrag news: snipped-for-privacy@e34g2000pro.googlegroups.com...

[...]

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:

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.

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

2:

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?

1. 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:

- 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:

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:

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

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

Elsevier/Newnes,

text -

Here is everything I used to compute the estimate for JCH's test. ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20SysID%20JCH.pdf One can see it didn't take a lot of effort. I tried n=3 first. It was easy to see the model required more poles. I tried n=5 next. The MSE is so small that it is close to the level of the noise. Minerr is a optimizing function that changes its parameters so MSE is minimized. The variable ERR is the returned MSE error that is used to provide the estimate of how good the fit it.

Peter Nachtwey

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

[...]

See PD5(PID)Z1Z2 control:

Page 1: State Observer Fig. 1 Page 2: Control without disturbance compensations Page 3: Control with disturbance compensations Page 4: Process Transfer Function F1(s) [=F3(s)], lower range

YhVvPBJULnanZ2dnUVZ snipped-for-privacy@comcast.com...

Why did you change the plant time constant from 0.123 to 0.006518124? Why don't you show your calculations. Anyone can generate lines. What happen if the feedback resolution is 0.0001?

Peter Nachtwey

Peter:

Why did you change the plant time constant from 0.123 to 0.006518124?

Jan:

It fitted to the bench mark defined earlier dicussed. Just to save some time rearranging the demo.

See for T=0.123

Page 1: State Observer Fig. 1 Page 2: Control without disturbance compensations Page 3: Control with disturbance compensations Page 4: Process Transfer Function F1(s) [=F3(s)], lower range

Peter:

What happen if the feedback resolution is 0.0001?

Jan:

T=0.0001? I would change to 0.1ms interpretation if 0.0001s was meant.

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

Not time, the feed back resolution. This means that the controller doesn't see perfect data. It sees data as if it came from a counter or AtoD converter and those counts must be scaled to provide the feed back data.

Peter Nachtwey

If you have an internal range 0.2 to 1 then 0.0001/1*100% = 0.01% resolution ist very good. If for some reason it were not then you can't make mathematics responsible for that.

The approximation error data chosen were +/- 1% max.

Identify this this one: ftp://ftp.deltacompsys.com/public/NG/plot1.zip The data is a CSV file format that can easily be imported into excel. The first column has units of seconds The second column is positions units in millimeters. The third column is the control output in volts. There are no more than 3 poles. There is no noise but the feed back is assumed to have a resolution of

10 microns. Have fun

Peter Nachtwey

"pnachtwey" schrieb im Newsbeitrag news: snipped-for-privacy@q5g2000prf.googlegroups.com...

Your data look like closed loop data (v2=volts changes).

I could have fun if I had just the process transfer function data like:

Number of data points = 20 would be ok. Just 10% step change of v2

t/seconds (e.g.) v1/millimeters (e.g.)

0 0 . . . . . . 0.48 .

It is open loop data. The system is also linear except for the small amount of quantizing that simulates a feed back resolution of 10 microns. .

You need to make your arrays bigger. You could use only every 20th or 25th point but that would not represent the control signal or the response very well. I doubt you can get a very good system identification if you do. Sampling must be done at small intervals in order to detect short time constants.

Exciting the system with step changes are not always permissible. From time to time we control systems in the 50 to 100 ton range and we have control much bigger loads. A step change can be very hard on machinery and excite the manager or owner more than the load. The system identification need to be able to handle an arbitrary control signal. This also makes it easy to identify systems with dead time.

Here is a hint. The answer is similar to the answer in the pdf file below. However, there is a 0.5 volt offset which causes the position to drift in the negative direction when the control output is 0. This simulates a valve being out of null or have an offset of half a volt. This is very realistic.

ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20Sysid2A2BV70%20T02.pdf

This is an example of system identification of a real hydraulic system. I have posted a link to this pdf before. The data is taken every millisecond for 1.5 second so there are 1500 points. Notice that the control output is also not a simple step change but has many small changes instead. I ploted the velocities instead of the positions because it is easier to see the results but feedback was in position units of inches..

Think of this problem as a puzzle. Have fun.

Peter Nachtwey

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

Sorry Peter, I don't like puzzles. I like solving problems like described below:

State of art high precision process identification and control tools where no puzzling is necessary.

Example: Input step v2=1, output differential equation (y=v1)

Testdata: e.g.

1y''''''+6y'''''+15y''''+20y'''+15y''+6y'+1y=1 F(s)=(1s+1)^6

Result:

1,000114y''''''+6,000568y'''''+15,00114y''''+20,00114y'''+15,00057y''+6,000113y'+1y=1,000063

F(s)=(1,000019s+1)^6

Error: 0,0019% after least square approximation of data. In practice 1% would be good enough.

But last not least, all tools must be easy to use! That means without instruction manual.

For new readers:

benchmark scheme and automatic PID-optimisation]

Newsbeitragnews:0KudnYexXZpKZ6ranZ2dnUVZ snipped-for-privacy@comcast.com...

???????

What if the input is not a step? Then your "state of high art precision process identification" isn't very useful.

This is not very realistic. For starters, you should be able to handle

1y''''''+6y'''''+15y''''+20y'''+15y''+6y'+1y=Gain*ControlOutput(t- DeadTime)+Offset

A gain is required. Hardly any system have a gain of 1 The control output function can be any given data. It can be sampled, a forcing function like a sine wave or the control output from a PID. The dead time is used to offset the time index into the forcing function. Systems tend to drift about in open loop mode so a offset is required. The pole locations can all be different. This means there are no simple binomial relationships between the constants.

In addition to the A0...A5 time constants, one must also determine the system gain, dead time and offset.

Tools must be useful too. If is up to the programmer ( you ) to make it easy.

Peter Nachtwey

I had never problems having a step in practice and then watch what was happening?

You can choose internally a gain of Kc = 1

Example:

Transmitter: Range 0...100°C

vT = 1/100°C

Kc = 1 (controller: 100%output/100%input -> P-band=100%)

Process: Range 0...10m^3

vP = 10m^3/1

All over chain gain v = vT*Kc*vP = 0.1m^3/°C

That means e.g. you need 0.1m^3 water for changing 1°C temperature.

See again Page 2, enlarged span lower figure (if reloading Fig.5)

is enough dead time that can be controlled.

Any other alternative is not available. That means you have no chance getting a better result.

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

The system is linear. It is open loop data. The control output is the default wave form we use to excite the system when auto tuning. Since we are often tuning hydraulic systems the process must be done extending and retracting because the surface area of the piston is different on the two sides so the gain is different in each direction.

I won't do that. Sometimes step changes aren't allowed. Imagine a

10% step change when moving 100 metric tons. That would excite the owners and managers and would test the welds and the foundations. The control input can be arbitrary becuse it often must be modified to suit the application.

In real life the people using the auto tuning don't know what the transfer function is. The auto tuning program must be able to try first order and second order under damped, critically damped and over damped etc. You must try all the models and select the best one.

In real systems one must be able to sample as fast as possible to identify those systems with very short time constants. Small motors often have time constants in the 3-5 millisecond range so fast sampling is required. I have said before we usually use about 500 to

2000 points. If 500 points are too much then you can use every 20th or 25th one but you can see that will not work well because you can't represent the control signal or the motion profile accurately using every 20th or 25th point. I think it is time to expand your arrays.

Have fun.

Peter Nachtwey

You use a wave form and I use a step, and find the process transfer function F1(s). The rest is mathematics like

Just for new readers:

benchmark scheme and automatic PID-optimisation]