Finding the deviation point from an initial linear trend with python - Stack Overflow

admin管理员组

文章数量:1123430

I am working with stress-strain curves which initially has a linear trend after which there the trend is variable from test to test. I am trying to find the point at which the data deviates from this initial linear trend - (pointed highlights by red arrow in image).

So far I have this code:

import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

df_input = pd.ExcelFile('data.xlsx')

sheet_names = df_input.sheet_names
dict_of_sheets = {}
for sheet in sheet_names:
    dict_of_sheets['data' + sheet] = pd.read_excel(df_input, sheet_name=sheet)

for name,df in dict_of_sheets.items():
    df['eff CP'] = df['CP'] - df['PP']
    df['axial stress'] = df['eff CP'] + df['Sd']

def fit_func(x, a, b):
    return a*x + b

fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(10, 8))

for i, (name, df) in enumerate(dict_of_sheets.items()):
    row = i // 2
    col = i % 2
    y = df['axial stress']
    x = df['Axial Strain']

    max_y = y.max()
    max_row = df.loc[df['axial stress'] == max_y]
    max_x = max_row['Axial Strain'].values[0]

    label_peak_og = 'Peak stress: {0:.1f} MPa'.format(max_y)

    short_df = df[x< 0.15]
    stress= short_df['axial stress']
    strain = short_df['Axial Strain']

    popt, pcov = curve_fit(fit_func, strain, stress)
    extended_strain = np.linspace(0, 0.2, 100)

    axs[row, col].plot(x, y,zorder=0, label='edited data')
    axs[row, col].set_title(name)
    axs[row, col].set_xlabel('Strain')
    axs[row, col].set_ylabel('Stress')
    axs[row, col].scatter(max_x, max_y, color = 'hotpink', marker='^', zorder=1, label=label_peak_og)
    axs[row, col].plot(extended_strain, fit_func(extended_strain, *popt), '--', label='Fitted function',zorder=3,)
    axs[row, col].legend()

plt.tight_layout()
plt.show() 

I attempted to calculate the residuals between the data and the fitted linear trend - but because the linear trend is based on a smaller data set I can't calculate the residuals with the original data as the arrays are different sizes.

I am working with stress-strain curves which initially has a linear trend after which there the trend is variable from test to test. I am trying to find the point at which the data deviates from this initial linear trend - (pointed highlights by red arrow in image).

So far I have this code:

import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

df_input = pd.ExcelFile('data.xlsx')

sheet_names = df_input.sheet_names
dict_of_sheets = {}
for sheet in sheet_names:
    dict_of_sheets['data' + sheet] = pd.read_excel(df_input, sheet_name=sheet)

for name,df in dict_of_sheets.items():
    df['eff CP'] = df['CP'] - df['PP']
    df['axial stress'] = df['eff CP'] + df['Sd']

def fit_func(x, a, b):
    return a*x + b

fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(10, 8))

for i, (name, df) in enumerate(dict_of_sheets.items()):
    row = i // 2
    col = i % 2
    y = df['axial stress']
    x = df['Axial Strain']

    max_y = y.max()
    max_row = df.loc[df['axial stress'] == max_y]
    max_x = max_row['Axial Strain'].values[0]

    label_peak_og = 'Peak stress: {0:.1f} MPa'.format(max_y)

    short_df = df[x< 0.15]
    stress= short_df['axial stress']
    strain = short_df['Axial Strain']

    popt, pcov = curve_fit(fit_func, strain, stress)
    extended_strain = np.linspace(0, 0.2, 100)

    axs[row, col].plot(x, y,zorder=0, label='edited data')
    axs[row, col].set_title(name)
    axs[row, col].set_xlabel('Strain')
    axs[row, col].set_ylabel('Stress')
    axs[row, col].scatter(max_x, max_y, color = 'hotpink', marker='^', zorder=1, label=label_peak_og)
    axs[row, col].plot(extended_strain, fit_func(extended_strain, *popt), '--', label='Fitted function',zorder=3,)
    axs[row, col].legend()

plt.tight_layout()
plt.show() 

I attempted to calculate the residuals between the data and the fitted linear trend - but because the linear trend is based on a smaller data set I can't calculate the residuals with the original data as the arrays are different sizes.

• python
• linear-regression

Share Improve this question Follow edited 13 hours ago travgeol asked 16 hours ago travgeoltravgeol 5944 bronze badges 7

• Any reason why you can't share your data? Stress-strain information is hardly subject to export controls. A couple of possibilities: (i) least-squares fit of straight line + parabola up to the maximum stress (which should be identifiable); (ii) smooth/filter your data, then look at where the derivative deviates from a constant. Are you looking at purely plastic materials or are you also considering brittle ones (which won't behave like that curve). – lastchance Commented 15 hours ago
• Hi @lastchance. I can't share the data as it is confidential but all the datasets generally shows a typical brittle response - so initially linearly elastic, then poroelastic behaviour with strain softening (though a few cases of strain hardening). In any case, here I am only really interested in the pre-yield behavior. – travgeol Commented 15 hours ago
• @lastchance I will have a go at finding how to implement your possible solution :) – travgeol Commented 15 hours ago
• Can you share some data "like" your data? Or even ... every second point with a multiplying factor? – lastchance Commented 15 hours ago
• Take first N points and build two regressions - linear and quadratic. Check R2 of each of the regression. Take difference of 'R2 linear' minus 'R2 quadratic' and find deviation point based on this difference value – Johnny Cheesecutter Commented 14 hours ago

 |  Show 2 more comments

1 Answer 1

Sorted by: Reset to default 0

Well, in the absence of your data, I "borrowed" some from a University website, https://efcms.engr.utk.edu/ef105-2021-01/grav/labs/optimization-1/stress-strain-example

I have cut the data off at the maximum stress. (What happens after that is probably time-dependent!) Then I put strain and stress as the first two columns of file data.csv

Then I fitted the data piecewise with a straight line y=Ex for x<s0 and y=ax²+bx+c for x>=s0. Here, parameters s0, E and a are to be found, whilst b and c are found from the condition that y and dy/dx are continuous at the join.

For this data it gave linear elastic limit s0=0.0202 and Young's modulus E=239984 (in whatever units they are using for their measurement).

The linear elastic limit isn't that critical; it could vary quite a lot. I guess you are more interested in the Young's modulus.

You could try other curves for your data if a parabola doesn't fit after the linear-elastic limit. Make sure that you cut off the data at maximum stress, though - what happens after that point is likely to depend on how you did your experiment. Irrelevant for brittle materials, which might actually break before they exhaust the linear regime.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

xdata, ydata = np.loadtxt( 'data.csv', skiprows=0, unpack=True, delimiter=',' )

def stressStrain( x, *params ):
    s0, E, a = params                  # x is strain, s0 is linear-elastic strain limit, E is elastic modulus (stress/strain)
    b = E - 2 * a * s0                 # ensures continuity and slope-continuity at elastic limit
    c = a * s0 ** 2                    #
    return np.piecewise( x, [ x < s0 ], [ lambda X: E * X, lambda X: a * X ** 2 + b * X + c ] )

p0 = ( 0.01, 1e6, 0 )
popt, pcov = curve_fit( stressStrain, xdata, ydata, p0=p0 )
s0, E, a = popt
print( 'Elastic limit,  s0: ', s0 )
print( 'Elastic modulus, E: ', E  )
print( popt )


# Plot data
xfit1 = np.linspace( 0 , s0       , 100 );   yfit1 = stressStrain( xfit1, *popt)
xfit2 = np.linspace( s0, xdata[-1], 100 );   yfit2 = stressStrain( xfit2, *popt)
plt.plot( xdata, ydata, 'k-', label='Data'           )
plt.plot( xfit1, yfit1, 'r-', label='Elastic region' )
plt.plot( xfit2, yfit2, 'b-', label='Plastic region' )
plt.xlabel("strain")
plt.ylabel("stress")
plt.legend()
plt.show()

Output:

Elastic limit,  s0:  0.020208648979847386
Elastic modulus, E:  239984.11710434384
[ 2.02086490e-02  2.39984117e+05 -2.33969193e+06]

本文标签： Finding the deviation point from an initial linear trend with pythonStack Overflow