【机器学习】python实现离散点拟合曲线并求曲率半径

目录

一、背景描述

二、问题描述

三、解决方案

---

一、背景描述

曲率半径是一种用来表征曲线上某处弯曲程度变化的量度，是一种灵敏度的表达形式，并且能够描述系统的平衡性状态。从数据驱动角度可知，数据变化幅度越大，曲率半径越小，系统平衡性越差；数据变化幅度越小，曲率半径越大，系统平衡性越好。 
当电网运行在稳定状态时，电网状态数据变化幅度较小，且位于合理区间。电网遭受扰动时，电网运行状态容易发生改变，电网响应的状态数据变化幅度较大。可见，系统平衡性与稳定性的变化，将伴随电网响应数据的波动。从数据驱动的角度可定义电网的平衡状态和不平衡状态：平衡状态下数据分布于合理的区间；不平衡状态下一定量数据分布于合理区间以外。多维度的数据可以用状态空间中主要的特征量及其变化来刻画系统的平衡性。曲率半径从数学意义上表征了系统的这种物理状态，考虑控制系统平衡性理论，本文将曲率半径作为系统平衡与不平衡的评价指标，当系统数据平衡性差时，对应数据曲率半径小；当系统数据平衡性好时，对应数据曲率半径大。综上可知，当系统数据曲率半径较小时，可判定为薄弱环节。

电网各元件之间有较强的耦合关系，一处扰动会引起其环节状态发生变化，将复电压数据作为电网薄弱环节统一分析方法的基础数据，复电压数据是电网运行过程的状态数据，其数据的波动性与平衡性，是电网运行状态的直观映射，复电压状态数据的研究对电网安全稳定分析有重要作用。本文采用节点复电压数据的曲率半径作为衡量其灵敏度变化的量度，即曲率半径作为节点复电压数据的变化特征。 

 

---

二、问题描述

需要根据已有的离散点、拟合出曲线函数，然后对每个点所在的地方求其曲率半径。

---

三、解决方案

代码：

import numpy as np
import matplotlib.pyplot as plt

#生成数据
x = np.linspace(0, np.pi, 10)
y = 5.0*np.cos(x)

#拟合曲线
fx = np.poly1d(np.polyfit(y,x,5))
dfx = fx.deriv()  # deriv()方法可得该格式函数的导函数
ddfx = dfx.deriv()

#求曲率半径
def calc_kappa(s):
    dx = dfx(s)
    ddx = ddfx(s)
    r = ((1 + (dx ** 2)) ** (3 / 2)) / ddx
    return r

snew = np.linspace(0, 10, 50)
k_buf = []
for i in snew:
    k = calc_kappa(i)
    k_buf.append(k)
k_buf = np.array(k_buf)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(snew[1:], k_buf[1:], s=1.0, color='red')
print(k_buf)
plt.show()

结果：

[-6.53193563e+15  1.33313971e+02  6.81808541e+01  4.72685489e+01
  3.75737277e+01  3.26003863e+01  3.03470185e+01  3.02572039e+01
  3.27199881e+01  3.98667234e+01  6.16589764e+01  2.78196094e+02
 -7.62336205e+01 -2.92568474e+01 -1.67384408e+01 -1.11393142e+01
 -8.06852782e+00 -6.19380404e+00 -4.97862784e+00 -4.16929732e+00
 -3.63384573e+00 -3.30058795e+00 -3.13163381e+00 -3.11074516e+00
 -3.23787238e+00 -3.52716870e+00 -4.00701257e+00 -4.72131964e+00
 -5.73179385e+00 -7.12096897e+00 -8.99600601e+00 -1.14932740e+01
 -1.47837723e+01 -1.90794686e+01 -2.46406346e+01 -3.17842696e+01
 -4.08937100e+01 -5.24295296e+01 -6.69418450e+01 -8.50841470e+01
 -1.07628790e+02 -1.35484278e+02 -1.69714487e+02 -2.11559997e+02
 -2.62461671e+02 -3.24086684e+02 -3.98357149e+02 -4.87481563e+02
 -5.93989251e+02 -7.20768017e+02]

注意： 

 这里不画第一个是因为第一个点的曲率半径绝对值极大：-6.53193563e+15。如果画的话，就会是这个样子：

代码：

import numpy as np
import matplotlib.pyplot as plt

#生成数据
x = np.linspace(0, np.pi, 10)
y = 5.0*np.cos(x)

#拟合曲线
fx = np.poly1d(np.polyfit(y,x,5))
dfx = fx.deriv()  # deriv()方法可得该格式函数的导函数
ddfx = dfx.deriv()

#求曲率半径
def calc_kappa(s):
    dx = dfx(s)
    ddx = ddfx(s)
    r = ((1 + (dx ** 2)) ** (3 / 2)) / ddx
    return r

snew = np.linspace(0, 10, 50)
k_buf = []
for i in snew:
    k = calc_kappa(i)
    k_buf.append(k)
k_buf = np.array(k_buf)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(snew, k_buf, s=1.0, color='red')
print(k_buf)
plt.show()

结果： 

 这样就看不出其他点的变化了。因此在处理时，不画第一个点。如下图所示：