APDL Pseudoinverse Least Square Fit

Curve fitting is a frequently used tool in engineering. I wished my linear algebra teacher taught me pseudoinverse. As used in the previous blog post script, it computes the least-square curve fit for linear equations. This can come in handy for fitting 2D or even 9D variables.

For a simple 2D straight line, the equation is \( y = m \cdot x + c\). Writing it matrix form, we have:

\(y = \begin{bmatrix} x & 1 \end{bmatrix} \cdot \begin{bmatrix} m\\c \end{bmatrix}\)

With multiple data points in space (e.g. \((x_1,y_1), (x_2,y_2),\dots\) ) the equation looks like this...

\( \begin{bmatrix}y_1\\y_2 \\y_3\\ \vdots \end{bmatrix} = \begin{bmatrix} x_1 & 1\\ x_2 & 1 \\ x_3 & 1 \\ \vdots & \vdots \end{bmatrix} \cdot \begin{bmatrix} m\\c \end{bmatrix} \)

This matrix looks similar to \(y = A \cdot b\) where our

\(A = \begin{bmatrix} x_1 & 1 \\ x_2 & 1 \\ x_3 & 1 \\ \vdots & \vdots \end{bmatrix} \) and \(b = \begin{bmatrix} m\\c \end{bmatrix} \)

If \(A\) is a square matrix, we can inverse \(A\) to get \(A^{-1} \cdot y = b\). However when the \(A\) matrix is not square, it is over constrained. There are more equations than unknowns like the problem at hand above. That's where the Moore-Penrose inverse does it's magic. From Wikipedia, the pseudoinverse is written as:

\(A^{+} = (A^{T} \cdot A)^{-1} \cdot A^{T} \)

where \(A^{+} \cdot A = I \)

Rewriting the linear equation to solve for \(b\):

\((A^{T} \cdot A)^{-1} \cdot A^{T} \cdot y = b\)

To program it in Ansys efficiently, the equation is rearranged as:

\((A^{T} \cdot A) \cdot b\ = A^{T} \cdot y \)

APDL script
Here's an APDL script that works through it.

!! Known answers to shoot for
!! y = mm*x+cc 
mm = 12
cc = 34

! Setup A Matrix and fake data with known answers
/prep7
nval = 20
*dim, Amat,, nval, 2
*dim, y,, nval, 1
*do,ct,1,nval
    Amat(ct,1) = ct ! x variable
    Amat(ct,2) = 1
    y(ct,1) = (mm*ct+cc)*(1+0.01*rand(-1,1)) ! adds random noise
*enddo

! Calculates Transpose
*dim, AmatT,, 2, nval
*mfun,AmatT(1,1),tran, Amat(1,1)

! Computes LHS Matrix
*dim, lhs,, 2,2
*moper, lhs(1,1), AmatT(1,1), mult, Amat(1,1) ! LHS = \((A^{T} \cdot A) \)

! Computes RHS Matrix
*dim, rhs,, nval,1
*moper, rhs(1,1), AmatT(1,1), mult, y(1,1) ! RHS = \(A^{T} \cdot y \)

! Solves matrix
*dim, coef,, 2,1
*moper, coef(1,1), lhs(1,1), solv, rhs(1,1) ! LHS*coef = RHS, solve for coef

! Prints out Gradient & Intercept
*stat, coef

Resources
Here's the text file of the above script: [Link]
Short lecture on it on YouTube: [Link]

Comments

  1. VillarealMarch 17, 2021 at 10:58 PM

    Thank you a lot, you solve my big problem

    ReplyDelete
  2. SimonJuly 8, 2023 at 11:59 AM

    Hi,

    Thank you so much for this code. But I'm facing an issue. I have extracted a sparse stiffness matrix of a large structure and for some calculation I need to find the pseudoinverse of the stiffness matrix. The code what you have provided does not work for sparse matrix. Is there any way I could do it?

    ReplyDelete
    Replies
    1. JasonJuly 8, 2023 at 12:54 PM

      The sparse matrix extracted are not in array format. You could do an inverse on the stiffness matrix while still in APDL math using *LSENGINE, *LSFACTOR and *LSBAC such as described in: https://www.ansystips.com/2017/10/export-stiffness-matrix-from-ansys.html.

      Delete
    2. SimonJuly 10, 2023 at 1:33 PM

      Hi Jason,

      Thank you for your reply. But if my stiffness matrix is singular, will it inverse using *LSENGINE, *LSFACTOR and *LSBAC?

      Delete
    3. JasonJuly 10, 2023 at 3:21 PM

      Unfortunately if your stiffness matrix is singular the inverse will not work. I don't think pseudoinverse will work for you either.

      Delete

Post a Comment

Popular posts from this blog

ANSYS APDL Syntax Highlighting editor

Image
Notepad++ with APDL User Defined Language The editor of my choice is Notepad++  with the available User Defined Language Files for APDL. You can install it without administrative privileges via the zip file. The best part of it is, it's FREE! After installing Notepad++, go to "Language>Define Your Language..." then "Import" the XML file downloaded from the above link. Remember to restart Notepad++ so that the language changes will take into effect. Opening up any *.inp or *.ans files should automatically switch highlighting to APDL. I made some minor edits. Here's my XML file: LINK . I also heard Sublime Text and  Ultraedit  has more advance features but they aren't (totally) free.
Read more

ANSYS User Defined Results

Image
There is an abundant of options in ANSYS classic when one wishes to post process results. ANSYS workbench default pull down menu post processing options are more limited but they can still be accessed via the User Defined Results. One way not commonly used but can come in handy is as follows: Zeroth: Under Analysis Settings, there is "Output Controls" where you can toggle to "Yes" what you would like to save before the solution starts. This is like OUTRES in APDL. Output Controls First: After solving the model, click on Solution in the tree to highlight it. Solution Second: Click on Worksheet in the toolbar. Worksheet Third: In the worksheet, you will see list of results that are saved. Right click on it to create the User Defined Results. Create User Defined Results So here we have it. You could of look up the different expressions in the help document but I find this method of accessing the results convenient.  Example: Aspec...
Read more

Export Stiffness Matrix from Ansys

Image
It is sometimes useful to extract the mass and stiffness matrix from Ansys.     *SMAT, MatK, D, IMPORT, FULL, file.full, STIFF       *PRINT, matk, matk, txt Exporting mass matrix would be similar:       *SMAT, MatM, D, import, full, file.full, MASS The above script uses APDL Math to get the job done. (Please see previous post for another example). The ordering of the matrix is unfortunately not concurrently exported. To verify the sequencing is as expected, we will work to replicate a truss example in the  Finite Element Trusses course notes by Bob Greenlee. Figure 1: Truss Problem Setup Model Creation Script to create model: /prep7 !! Creates Model to reflect course notes ! Properties et ,1,1  mp , ex, 1, 29.5e6 r , 1, 1 ! Geometry n ,1 $  n ,2, 40 $  n ,3, 40, 30 $  n ,4, 0, 30 e ,1,2 $  e ,2,3 $  e ,1,3 $  e ,3,4 ! Boundary Conditions d ,1,ux,0 $  d ,1,uy,0 d ,2,uy...
Read more