PD Model: Curve fitting with Python

With a set of data, we will see how we can fit the data the phenomenological model by using \(A\) and \(b\) parameters with Python. A sample PD data can be downloaded on Github. The curve fitting process in Python is fairly easy. It uses a non-linear least square method to fit a function. We will need to import the Scipy libraries to perform the fitting. The code snippet is provided as below:

Lines #1-4: Read the .csv file. We see that there's a total of 251 lines of data with each line contains permanent strain in percent, number of load cycle, deviatoric stress, and shear stress ratio. Obviously, in this exercise, we are only interested in the accumulation of permanent strain due to load cycles. Therefore, we are going to ignore the effects of stress state for now.

	P_Strain_%	Load_Cycle_N	dS_psi	SSR
0	0.17	1	16.03	0.25
1	0.32	41	16.19	0.25
2	0.35	81	16.11	0.25
3	0.37	121	16.27	0.25
4	0.38	161	16.11	0.25
...	...	...	...	...
246	0.52	9841	16.19	0.25
247	0.52	9881	16.11	0.25
248	0.52	9921	16.19	0.25
249	0.52	9961	16.19	0.25
250	0.52	9996	16.03	0.25

Line #7: explores the correlation between each feature. This is a quick way to see any correlation between each feature or variable. For example, we see that the strongest relationship exists between P_Strain_% and Load_Cycle_N, i.e. 0.792972. Another strong positive relationship also exists between dS_psi and SSR. This is expected because deviatoric stress (dS_psi) is directly computed from shear stress ratio derived from the Mohr-Coulomb envelope. Again, we will come back to the effects of stress state later.  

	P_Strain_%	Load_Cycle_N	dS_psi	SSR
P_Strain_%	1.000000	0.792972	-0.096899	-0.036616
Load_Cycle_N	0.792972	1.000000	-0.079329	-0.040335
dS_psi	-0.096899	-0.079329	1.000000	0.764572
SSR	-0.036616	-0.040335	0.764572	1.000000

Lines #9-29: This is where we define the phenomenological model and run the curve fitting method. Note that we call "curve_fit" from the scipy.optimize library, then fit the "phe_model" with the xdata and ydata, which respectively, are number of load cycle and permanent strain. Note that "curve_fit" includes several other useful parameters, including uncertainty (sigma), initial guess (p0), and bounds etc. You do not need to input these parameters if you do not know any. The breakdown of the next few lines are as below:  

• Line #22 returns an array of the model parameters \(A=0.25561268\) and \(b=0.0786587\). The \(b\) value is reasonable as discussed in the previous post. 
• Lines #23 &24 provides a covariance matrix with respect to \(A\) and \(b\); all of which are near zero. 
• Lines #27-29 is a further step to find out the \(R^2\) value of the fitting. In this case, \(R^2 = 0.965\). This high \(R^2\) value should not be a surprise for a power relationship.

Now, let's plug back \(A\) and \(b\) to the model and plot the graph. An easy way is to create an empty array, p_strain, using Numpy and parse in \(A\) and \(b\) to the "phe_model" function. The computed permanent strain are stored in the p_strain array. Following that, we plot both original data and model results against the number of load cycle.

Read more about scipy.optimize.curve_fit

Download sample file from Github