Polynomial fitting on Mitsubishi FX5u

[luca72]

Yes there is: Writing code to do Gaussian Elimination (actually Gauss-Jordan; see here; there is a flowchart that may be helpful for converting this to Ladder Logic) of a 4x4* array is probably the easiest way to replicate, on the PLC, the functionality of the Python numpy.polyfit routine to calculate the third-order polynomial coefficients to model a least-squares fit to the 5 to 10 [position, radius] pairs. It will be much easier with Structured Text, if that is available, but it can still be done with Ladder Logic. There are other, faster ways to factor the array, but they are more complex and would be harder to code and debug. With Gaussian Elimination, there are only a few operations: scaling a row; adding one (scaled) row to another; swapping rows, It might not be necessary to swap rows. The only problems will be an ill-conditioned array and numerical instability.

The time we lost waiting for the answer is small; the time you will spend, or that you pay someone to spend, coding that solution will not be small.

* When I wrote "3x3" earlier that was a mistake; I was forgetting the zeroth-order term.

thanks for your reply what i what to do is replicate in st but for me is not easy and for have the same result i need to do all the stuf that polyfit do
thanks

 

[drbitboy]

thanks for your reply what i what to do is replicate in st but for me is not easy and for have the same result i need to do all the stuf that polyfit do
thanks

You do not need to do all the stuff that polyfit does: polyfit calls numpy.linalg.lstsq, and .lstsq uses singular value decomposition for efficiency because people are calling it with large matrices.

You have a 4x3 matrix (4x5 when augmented), and you want updated polynomial coefficients after several manual measurements, so (1) your matrix is small, and (2) you can wait for the answer i.e. you don't need the calculations to complete in one scan cycle (although they probably would).

So I looked online for "pascal gauss jordan" and I found this; that .7z (zip) file contains code that will be very similar to the ST you need. You still need to build the array.

 

[Peter Nachtwey]

Guass-Jordan, really? That will be error prone and fun to support. What happens when the OP changes the number of points to 10? He said there could be 5 to 10 points above. Think about the OP has to work with. However, that is the main problem. Also, it seems the OP is obsessed with finding a least squares fit rather than a cubic that goes through all the points. Python has a CubicSpline function that is easy to use but the OP keeps using a least squares function.

 

[drbitboy]

Phhht.

Third-order Gauss-Jordan is maybe several dozen ST statements; not error prone at all (except for ill-conditioned matrices, and the code can check for that) and ask for another set of matrices. The only issue might be if the FX does not have double-precision floating-point, but that is only maybe.

Once it's debugged and working there is nothing to support; it just works (or fails due to ill-conditioned matrix).

Cubic fit that does not have to pass through any of the measurements is what they want: Python's polyfit gave a satisfactory answer; duplicating it is time consuming but not complex.

The number of points is irrelevant as long as they know the upper limit beforehand; the solution approach is the same.

Spline or Lagrange will be better in my mind, but not be any less complex: whatever gain in simplicity comes with the algorithm is lost in the code that selects which samples are used in the evaluation.

 

[drbitboy]

Despite @Peter Nachtwey's very reasonable concerns, Gauss-Jordan works well as I expected, at least for this limited application.

See the image below for the results of a hard-coded PLC-based ST* solution that calculates cubic coefficients compared to the Python numpy.polyfit result (Siemens on the left; Python in the terminal window with the black background on the right).

The Python code uses double-precision REALs, while the PLC code** uses single-precision REALs, so there are some differences, but the coefficients agree to within a percent and the calculated evaluations agree to ~0.1; those differences can be attributed to the differences in floating-point precision, so the Python and PLC versions can be said to produce essentially identical results. I also ran the PLC code with double-precision LReals and saw agreement out to several more decimal places.

The point is that the matrix solver used (Gauss-Jordan; SVD; etc.) is unimportant; the generated polynomial is the same for all solver because the least-squares solution is unique. If double-precision REALs are not available in the FX, I think there are some scaling techniques that could be applied to reduce the roundoff/instability errors.

* SCL in Siemens TIA ≡ Structured Text
** Siemens TIA/S7-1200

The Python code is simple, as described by the OP:

% cat test_polyfit.py
import numpy

def test_cubicfit(xsamples, ysamples):
  c = numpy.polyfit(xsamples, ysamples, 3)
  return c

if "__main__" == __name__:
  import pprint
  ysamples = [10250, 9870, 8536, 7653, 5560]
  xsamples = [350, 460, 530, 680, 900]
  coeffs = test_cubicfit(xsamples,ysamples)

  print('\n=======XY pairs')
  for pair in zip(xsamples,ysamples): print(pair)

  print('\n=======Cubic coefficients:')
  print(coeffs[::-1].reshape((-1,1,)))

  print('\n=======Calculated Y values:')
  print(numpy.polyval(coeffs,xsamples).reshape((-1,1,)))
%

​

 

[Peter Nachtwey]

If you want least squares then the Savitsky-Golay, least squares, algorithm I posted many posts ago is the way to go. There are other methods that will provide the same results. I know that I used the S-G algorithm for autotuning in the RMC100 25 years ago. The points don't need to be equally spaced but they are in my example. Use wxMaxima to enter the array SYMBOLICALLY! It will generate the equation needed for that many points. This will need to be done once for each number of points the OP has. Our RMCTools also used the S-G algorithm to smooth graphical data.
look at this again
Savitzky-Golay
%O4 has the coefficients for the cubic equation y=a+bx+cx^2+dx^3
Notice that these coefficients also provide the first, second and third derivatives at any point. I can use these for feed forwards.
There is no indexing. It may be tedious to enter the equations for the 4 coefficients but the is no indexing and the computation time is deterministic.
This would need to be done for all the different cases but once entered he is done.
Indexing on a PLC is not easy or fast.

I would still use Lagrange interpolation. It isn't as difficult as it looks. This is my trick to avoid array indexing. Suppose there are 10 points and you need to interpolate between points 3 and 4. I would use the PLC's copy command to copy command to copy points 2,3,4,5 to a static area that is HARD CODED to compute the interpolated value. When the index moves to between points 4 an 5 I would copy the points 3,4,5 and 6 to the same static area and run the hard coded code. By copying to a static area you avoid lots of index. Likewise for what you want to do. Instead of indexing every variable. Copy the row and column that need to be multiplied to a static area and run the hard coded code. This reduces the indexing significantly and is faster.

If I were engineering manager I would ask how many machines all this programming can be amortized over. If it is just one or two I would say buy a motion controller with cubic splines. I would want to know why a least squares fit is better than a accurate fit. Lagrange interpolation with my static memory trick is so simple. On top of all of that I would want to know if the rotation is cyclic or does it just go between two angles

This is another fail thread. The main problem is using an under powered PLC.

@drbitboy
Edit, how do you do that in a PLC? PLCs can't use numpy or scipy. This is the same problem I had with writing firmware for the motion controller. I had to write all the code without using a library.

 

[busarider29]

This doesn't much help the OP (Peter and Brian have provided solutions, like usual), however this could be useful for others in the future, or just FYI. There is a matrix math library available (free) for download for TwinCAT. Someone was ambitious enough to do the hard work in creating it. I've only just downloaded and installed it, and have not tested out any of the FBs. However, I have an application where this is probably what I may need to use, and the OP probably could have used something like this too, although he/she is not using TC.

GitHub - BurksEngineering/TcMatrix: Matrix Math Library for Beckhoff TwinCAT 3 PLC Environment
Releases · tcunit/TcUnit

 

[drbitboy]

@drbitboy
Edit, how do you do that in a PLC? PLCs can't use numpy or scipy. This is the same problem I had with writing firmware for the motion controller. I had to write all the code without using a library.

I did the same: I hard-coded it.

Gauss-Jordan is O(N³), and N=4 for cubic, so it's ~64 (=4³) operations; hardly the labor of Hercules.

 

[drbitboy]

This doesn't much help the OP (Peter and Brian have provided solutions, like usual), ...

And now you have provided a solution too, only OP can't use yours either . I'm waiting for OP to offer me some Lira .

Welcome to the club; we're having jackets made up.

 

[Peter Nachtwey]

And now you have provided a solution too, only OP can't use yours either . I'm waiting for OP to offer me some Lira .

There have been too many posts lately where the OP is clueless as to what he really needs.
There was a recent thread about a sawmill application with a supposed 50 ft cylinder. They don't exist and if they did, they wouldn't work well.
I am in old fuddy duddy mode. No non-sense. If the OP has an engineering manager I would fire him. Next in line is the OP for wasting our time by not being up front with his requirements.
Here is a video from 2006 and one of our customers in Canada. This is a second-generation motion controller.
A new cubic spline is being downloaded while the previous cubic spline is being executed. Basically, the cubic splines are updated on-the-fly.

https://deltamotion.com/peter/Videos/JAN-04%20VSS_0001.mp4

Perhaps the OP should contact Tinine. Remember him. I made pipe bending controls with a proper chip too. He gave me lots of sh!t. However, I doubt he could get his propeller chip to do the required math out of the box.

BTW, I see this all the time on engineering forums. People want to use Matlab to get answers but are totally stumped when it comes to embedding algorithms in firmware code.

 

[busarider29]

BTW, I see this all the time on engineering forums. People want to use Matlab to get answers but are totally stumped when it comes to embedding algorithms in firmware code.

I agree, but then again I think it's about speed to development and deployment. I've got enough math skills and certainly enough programming background, a few university STEM degrees (including AI/data science), years of experience in "controls", yet all of that was still not enough to land the position I recently interviewed for. I'm almost certain it was because I didn't have MATLAB experience, as that is what they were using to "code" algorithms and flash to their microcontrollers. But I digress.

 

[Peter Nachtwey]

I agree, but then again I think it's about speed to development and deployment.
I've got enough math skills and certainly enough programming background, a few university STEM degrees (including AI/data science), years of experience in "controls", yet all of that was still not enough to land the position I recently interviewed for. I'm almost certain it was because I didn't have MATLAB experience, as that is what they were using to "code" algorithms and flash to their microcontrollers. But I digress.

That is a recent ability of Matlab. Python and Mathematica can do it too. You can get your app to run in Matlab and then convert it to C. I don't know how well it works and if it allows you to include Matlab's libraries compiled in C. The same can be done for python. I also write code in python that generates ST code for our motion controller. I must only replace the = with := but that is easily done with an editor.

Scilab is a lot like Matlab. If you can't afford Matlab then try Scilab. There are a few differences but if you learn Scilab then Matlab will be easy

There is always something to learn.

 

[drbitboy]

thanks for your reply what i what to do is replicate in st but for me is not easy and for have the same result i need to do all the stuf that polyfit do
thanks

Look for "gauss jordan reduced row echelon form" on youtube. that is what you need to do.

 

[busarider29]

That is a recent ability of Matlab. Python and Mathematica can do it too. You can get your app to run in Matlab and then convert it to C. I don't know how well it works and if it allows you to include Matlab's libraries compiled in C. The same can be done for python. I also write code in python that generates ST code for our motion controller. I must only replace the = with := but that is easily done with an editor.

Scilab is a lot like Matlab. If you can't afford Matlab then try Scilab. There are a few differences but if you learn Scilab then Matlab will be easy

It can be done for ST code too - build a MATLAB Simulink model and convert it to ST in Matlab. I want to get there but I'm not fluent in Matlab.....yet. Company purchased me a Matlab (with Simulink) license, and I bought a book.

I agree, there's always something to learn. I finally start getting good with Python (ML applications), and now it looks like Rust will likely be the future for that.

 

[Peter Nachtwey]

Look for "gauss jordan reduced row echelon form" on youtube. that is what you need to do.

I haven't used those techniques since college. Now I have faster algorithms but I wouldn't want to implement them in ladder. I wouldn't even want to implement Gauss-Jordan in ladder.
I can't see where you have provided a solution that runs exclusively on a PLC.
That is a lot of matrix manipulation if the OP has 10 points or for a general solution for 5 to 10 points.

Pipe bending should be simple. Remember Tinine? He claimed to sell hundred of pipe bending machines a year. He was the one that talked about the glory of the propeller chips and how the RMC's were trash.
Tinine does pipe bending in the UK. I bet he could find an answer on a propeller chip.

We still don't know why a least squares fit is required instead of a cubic spline that goes through all the points ( knots ).
Anyway, I would have had this done long ago. I could also make it so the points can be updated on-the-fly.
Radius was mentioned. That means the application is rotational. We don't know if it is cyclic "periodic" or just moves between two angles. If the application is cyclic then it is easier yet.

It is time to pop some pop corn.