Python数据处理：应力-应变曲线插值分析及潜在错误处理

import pandas as pd
import numpy as np
from scipy.interpolate import interp1d

M0 = pd.read_csv('E:\shujushuchu\shuchu1.csv')  #导入四个参数
M1 = pd.read_csv('E:\shujushuchu\shuru1.csv')   #导入abaqus生成的应力应变曲线的点
B0 = M0.values  
B1 = M1.values  #将导入到参数转换为矩阵
print(np.size(B1, 0) / 2)  #0表示行数，1表示列数，并赋值给a1
print(B1[1, 1])  #因为py语言里下标从0开始，检查数据的位置是否正确
print(B1.shape)  #根据这个值来修改B3、B4的矩阵大小,将后者赋值给z
a1 = 1
z = 160
B3 = np.zeros((a1, z))  #储存曲线点的矩阵
B4 = np.zeros((a1, z))  #储存预测参数的矩阵
                          #这两个矩阵根据前一步的计算结果修改这两个矩阵的大小，a1做为行元素的大小，而列元素大小要大于B1的shape
a2 = 0  #用于循环填充B3

for i in range(1, a1 + 1):  # a1值加1放在这儿，总循环，以行为单位
    a3 = np.sum(~np.isnan(B1[i, :])) #表格到进来后存在nan值，此处去除nan值，并且获得该行的长度，之后用到a3值错位行的长度，才能够实现插值（存在0会报错）
    B5 = np.zeros((2, a3))  #B5矩阵用于两两储存对应的应力应变点，便于之后的循环
    for i1 in range(1, a3 + 1):
        B5[0, i1 - 1] = B1[2 * i - 2, i1 - 1]
        B5[1, i1 - 1] = B1[2 * i - 1, i1 - 1]   #将应力应变数据储存为两行数据

    B2 = np.zeros((2, a3))  #创建一个空矩阵（该矩阵将用于储存应力应变曲线塑性阶段的点）
    C = B5[1, 2] / B5[0, 2]  #计算弹性模量
    a5 = np.size(B5, 1)  #计算B5矩阵有多少行
    a4 = 2   #用于之后循环跳出
    B7 = np.zeros((2, 3))   #创建空的B7矩阵便于储存条件屈服点（从弹性阶段进入塑性阶段）附近的三个点（用于寻找条件屈服点）
    for ii in range(2, a5):  

        if B5[1, ii - 1] <= C * (B5[0, ii - 1] - 0.002):  #寻找生成的数据中从弹性阶段进入塑性阶段瞬间的后一个点
            a4 = a4 + 1
            if a4 == 3:
                a6 = ii
                B7[0, 2] = B5[0, ii - 1]
                B7[1, 2] = B5[1, ii - 1]
            B2[0, a4 - 1] = B5[0, ii - 1]
            B2[1, a4 - 1] = B5[1, ii - 1]  #这两行循环储存条件屈服点后的实验数据于B2矩阵中


    B7[0, 0] = B5[0, a6 - 3]
    B7[0, 1] = B5[0, a6 - 2]
    B7[1, 0] = B5[1, a6 - 3]
    B7[1, 1] = B5[1, a6 - 2]    #该部分在B7中储存了条件屈服那一瞬间的后一个点和前两个点（用于后面的插值计算）
    B2[0, 1] = B5[0, a6 - 2]
    B2[1, 1] = B5[1, a6 - 2]    #该部分在B2中储存了条件屈服那一瞬间的后一个点
    a11 = (B7[0, 2] - B7[0, 1]) / 500   #将储存的三个点长度细分单位长度，用于插值计算获得中间更多的点
    f1 = interp1d(B7[0, :], B7[1, :], kind='linear')  #定义插值计算的数组、范围以及插值类型
    for i2 in np.arange(B7[0, 1], B7[0, 2], a11):  #利用for循环将数据细分并一一进行插值计算

        a12 = f1(i2)
        a13 = a12
        if a12 <= C * (i2 - 0.002):   #寻找条件屈服点的准确大概值
            B2[0, 1] = i2 - a11         #将条件屈服点应变储存进B2
            B2[1, 1] = a13 - C * a11      #将条件屈服点应力储存进B2
            break

    B6 = np.zeros((2, 10000))       #创建空白矩阵，用于插值计算条件屈服点后的点
    a7 = (np.max(B2[0, :]) - B2[0, 1]) / 500    #将塑性阶段细分单位长度用于插值计算
    a10 = 0    #用于循环递增储存插值后的细分的塑性阶段的点
    B2 = B2[:, np.nonzero(np.sum(B2, axis=0))[0]]    #去除矩阵末尾没储存数据的0值，保证插值计算的正常进行
    f2 = interp1d(B2[0, :], B2[1, :], kind='cubic')    #定义插值计算的数组、范围以及插值类型
    for i3 in np.arange(B2[0, 0], np.max(B2[0, :]) + a7, a7):    #利用for循环将数据细分并一一进行插值计算
        a20 = i3 + a7
        a22 = np.max(B2[0, :])
        if a20 > a22:
            a20 = a22
        a9 = f2(a20)
        a23 = i3
        if a23 > a22:
            a23 = a22
        a8 = f2(a23)       #为了寻找塑性失稳点，差值范围之前设置只到了原始数据的最大值处，插值时可能会超过这个范围而不能正常进行，所以将超过的部分等最大点处的数据
        derivative = (a9 - a8) / a7   #计算曲线的斜率
        if derivative - a9 >= 0:    #在斜率减缓到0之前，将所有的插值点都循环储存到B6矩阵里
            a10 += 1
            B6[0, a10 - 1] = i3
            B6[1, a10 - 1] = a8


    
    B6 = B6[:, np.nonzero(np.sum(B6, axis=0))[0]]     #去除矩阵末尾没储存数据的0值，保证插值计算的正常进行
    df = pd.DataFrame(B6)
    writer = pd.ExcelWriter('output3.xlsx', engine='openpyxl')
    df.to_excel(writer, sheet_name='Sheet1', index=False)
    writer.book.save('output3.xlsx')

    
    if np.max(B6[0, :]) > 0.47:      #将未能插值计算到应力极限点的数据去除掉（不能进行接下来的循环的话将不会储存进最终的矩阵里）
        continue
    else:
        a2 += 1
        for iiii in range(1, 5):
            B4[a2 - 1, iiii - 1] = B0[i - 1, iiii - 1]    #将参数储存进B4矩阵中

        B3[a2 - 1, 9] = np.min(B6[0, :])
        B3[a2 - 1, 10] = np.max(B6[0, :])    #将塑性阶段应变的最大值与最小值储存到B3矩阵中的第9列与第10列
        B6[0, :] = (B6[0, :] - B6[0, 0]) / (max(B6[0, :]) - B6[0, 0])   #将B6中的数据做归一化处理便于之后的插值计算（用于插值时取对应比例处的值）
        
                
        f = interp1d(B6[0, :], B6[1, :], kind='cubic')    #定义插值计算的数组、范围以及插值类型
        xx = np.array([0, 0.03, 0.09, 0.18, 0.3, 0.45, 0.63, 0.84, 1.0])    
yy = f(xx)                                              #取得xx比例处的值

        for iii in range(1, 10):
            B3[a2 - 1, iii - 1] = yy[iii - 1]  #将按比例选取的9个值储存进B3矩阵当中
            
#完成数据的处理

B3 = B3[:, np.nonzero(np.sum(B3, axis=0))[0]]
B4 = B4[:, np.nonzero(np.sum(B4, axis=0))[0]]  #去除矩阵中的多余0的值便于储存数据
df = pd.DataFrame(B3)
writer = pd.ExcelWriter('output1.xlsx', engine='openpyxl')
df.to_excel(writer, sheet_name='Sheet1', index=False)
writer.book.save('output1.xlsx')
df = pd.DataFrame(B4)
writer = pd.ExcelWriter('output2.xlsx', engine='openpyxl')
df.to_excel(writer, sheet_name='Sheet1', index=False)
writer.book.save('output2.xlsx')                          #将数据生成为表格

上述代码为数据处理代码，一般情况下能正确的运行处理出正确的结果，但是处理很多数据的时候，会偶尔出现某组数据处理错误的情况，这是什么原因呀，是因为这段程序里面有缺陷吗，具体出现的处理错误的情况是插值寻找寻找塑性失稳点数值出现错误或者是f = interp1d(B6[0, :], B6[1, :], kind='cubic')    #定义插值计算的数组、范围以及插值类型
        xx = np.array([0, 0.03, 0.09, 0.18, 0.3, 0.45, 0.63, 0.84, 1.0]) 这段程序出现负数的情况，正确的情况应该是从前往后依次递增的数据，请中文解释一下，因为大部分情况都能正常运行，最主要数据不可能有问题，因为单独拿出来处理的时候就会正确得到结果，但是很多组数据一起算的时候就可能出现某些组处理出错的情况，所以代码是什么问题

根据您提供的代码，有几个可能导致处理错误的原因：

1. **数据缺失：** 在处理数据之前，您使用了`np.isnan()`函数来检查NaN值，并根据结果计算行的长度。如果数据中存在NaN值，可能会导致计算行长度出现错误。建议在处理数据之前，先对数据进行清洗，确保没有NaN值。

2. **数据索引错误：** 在代码中，您使用了`B1[2 * i - 2, i1 - 1]`和`B1[2 * i - 1, i1 - 1]`来获取应力应变数据的点。请确保索引值是正确的，以避免获取错误的数据点。

3. **参数矩阵大小错误：** 根据您的代码，参数矩阵B4的大小是根据参数矩阵B0的大小来确定的。请确保两个矩阵的大小匹配，否则可能会导致数据存储错误。

4. **插值计算错误：** 在使用`interp1d`函数进行插值计算时，可能会出现插值点超出原始数据范围的情况。您可以通过检查插值点的范围是否在原始数据范围内来解决这个问题。

5. **塑性失稳点识别错误：** 代码中使用斜率的变化来识别塑性失稳点，可能会受到数据噪声的影响。建议对数据进行平滑处理，或者使用更稳健的识别方法。

6. **插值函数类型错误：** 您使用了`kind='cubic'`的立方插值方法，这可能会导致插值结果出现波动，尤其是在数据点稀疏的情况下。建议尝试使用线性插值或其他更平滑的插值方法。

综上所述，建议您先检查数据是否完整，然后仔细检查索引和矩阵大小是否正确。另外，您还可以添加一些错误处理的机制，例如使用try-except语句来捕获异常并处理错误的情况。