使用python进行模糊控制模拟仿真
输入量为和
，单输出
构建一个模拟环境EPmodel.py，一个输入一个输出
import numpy as np
import matplotlib.pyplot as plt
import math
import randomPa = 101class EPmodel:def __init__(self):self.a = 973self.n = 0.955self.eta = 0.775self.diameter = 6.280 self.A = 1/4 * math.pi * self.diameter * self.diameterself.L = 900 self.Vc = self.A * self.L * 0.599self.As = 0.953 self.h = 960 def Pressure(self, p, v, n):K1 = p / PaK2 = self.a * math.pow(K1, self.n) / self.VcK3 = self.A * vK4 = (self.As * self.h * self.eta * n) / (2 * math.pi)K5 = K3 - K4deltaP = K2 * K5p_ = p + deltaPp__ = np.random.normal(loc=p_, scale=2.0, size=1)return p__[0]if __name__ == '__main__':#测试模拟环境epm = EPmodel()p = 120v = 40n = 15v = v / 60n = n / 60pList = []for i in range(10):p_ = epm.Pressure(p, v, n)p = p_print(p)pList.append(p)plt.plot(pList)plt.show()
编写模糊控制器controller.py
import numpy as npclass FuzzyController:def __init__(self, n0):self.outputLast = n0self.outputMin = 0self.outputMax = 16def TensorProduct(self, a, b):## a是列向量(表示有多少行) b是行向量(表示有多少列)rows = len(a)cols = len(b)c = np.ones((rows, cols))ChooseSmall = lambda x, y: y if x > y else xfor i in range(rows):for j in range(cols):c[i, j] = ChooseSmall(a[i], b[j])return cdef ReshapeToCol(self, a):rows = a.shape[0]cols = a.shape[1]length = rows * colsb = np.ones(length)for i in range(rows):r = i*rowsc = i*rows+colsb[r:c] = a[i, :]return b[..., np.newaxis]def ReshapeToRow(self, a):rows = a.shape[0]cols = a.shape[1]length = rows * colsb = np.ones(length)for i in range(rows):r = i * rowsc = i * rows + colsb[r:c] = a[i, :]return b[np.newaxis, ...]def MatricUnion(self, a, b):## 矩阵 a b 取并集c = np.ones((a.shape[0], a.shape[1]))for i in range(a.shape[0]):for j in range(a.shape[1]):if a[i, j] > b[i, j]:c[i, j] = a[i, j]else:c[i, j] = b[i, j]return cdef VectorIntersection(self, a, b):##两个向量取交集（取小数）ChooseSmall = lambda x, y: y if x > y else xc = np.ones(a.shape[0])for i in range(c.shape[0]):c[i]= ChooseSmall(a[i], b[i])return cdef Compose(self, a, b):##矩阵合成运算i = a.shape[0]j = b.shape[1]c = np.ones((i, j))for ii in range(i):for jj in range(j):c[ii, jj] = self.VectorIntersection(a[ii, :], b[:, jj]).max()return cdef Output(self, Err, ErrK):## Err是当前值-目标值 x - xd# -40, -20, 0, 20, 40eTable = [[1, 0, 0, 0, 0],[0, 1, 0, 0, 0],[0, 0, 1, 0, 0],[0, 0, 0, 1, 0],[0, 0, 0, 0, 1]]eTable = np.array(eTable)# -40, -20, 0, 20, 40ecTable = [[1, 0, 0, 0, 0],[0, 1, 0, 0, 0],[0, 0, 1, 0, 0],[0, 0, 0, 1, 0],[0, 0, 0, 0, 1]]ecTable = np.array(ecTable)# -4, -2, 0, 2, 4uTable = [[1, 0, 0, 0, 0],[0, 1, 0, 0, 0],[0, 0, 1, 0, 0],[0, 0, 0, 1, 0],[0, 0, 0, 0, 1]]uTable = np.array(uTable)R = np.zeros((25, 5))## 控制规则a1 = eTable[0, :]a2 = eTable[1, :]a3 = eTable[2, :]a4 = eTable[3, :]a5 = eTable[4, :]b1 = ecTable[0, :]b2 = ecTable[1, :]b3 = ecTable[2, :]b4 = ecTable[3, :]b5 = ecTable[4, :]ab11 = self.ReshapeToCol(self.TensorProduct(a1, b1))ab12 = self.ReshapeToCol(self.TensorProduct(a1, b2))ab13 = self.ReshapeToCol(self.TensorProduct(a1, b3))ab14 = self.ReshapeToCol(self.TensorProduct(a1, b4))ab15 = self.ReshapeToCol(self.TensorProduct(a1, b5))ab21 = self.ReshapeToCol(self.TensorProduct(a2, b1))ab22 = self.ReshapeToCol(self.TensorProduct(a2, b2))ab23 = self.ReshapeToCol(self.TensorProduct(a2, b3))ab24 = self.ReshapeToCol(self.TensorProduct(a2, b4))ab25 = self.ReshapeToCol(self.TensorProduct(a2, b5))ab31 = self.ReshapeToCol(self.TensorProduct(a3, b1))ab32 = self.ReshapeToCol(self.TensorProduct(a3, b2))ab33 = self.ReshapeToCol(self.TensorProduct(a3, b3))ab34 = self.ReshapeToCol(self.TensorProduct(a3, b4))ab35 = self.ReshapeToCol(self.TensorProduct(a3, b5))ab41 = self.ReshapeToCol(self.TensorProduct(a4, b1))ab42 = self.ReshapeToCol(self.TensorProduct(a4, b2))ab43 = self.ReshapeToCol(self.TensorProduct(a4, b3))ab44 = self.ReshapeToCol(self.TensorProduct(a4, b4))ab45 = self.ReshapeToCol(self.TensorProduct(a4, b5))ab51 = self.ReshapeToCol(self.TensorProduct(a5, b1))ab52 = self.ReshapeToCol(self.TensorProduct(a5, b2))ab53 = self.ReshapeToCol(self.TensorProduct(a5, b3))ab54 = self.ReshapeToCol(self.TensorProduct(a5, b4))ab55 = self.ReshapeToCol(self.TensorProduct(a5, b5))r11 = self.TensorProduct(ab11, uTable[0, :])r12 = self.TensorProduct(ab12, uTable[0, :])r13 = self.TensorProduct(ab13, uTable[0, :])r14 = self.TensorProduct(ab14, uTable[1, :])r15 = self.TensorProduct(ab15, uTable[3, :])r21 = self.TensorProduct(ab21, uTable[0, :])r22 = self.TensorProduct(ab22, uTable[0, :])r23 = self.TensorProduct(ab23, uTable[1, :])r24 = self.TensorProduct(ab24, uTable[2, :])r25 = self.TensorProduct(ab25, uTable[3, :])r31 = self.TensorProduct(ab31, uTable[0, :])r32 = self.TensorProduct(ab32, uTable[1, :])r33 = self.TensorProduct(ab33, uTable[2, :])r34 = self.TensorProduct(ab34, uTable[3, :])r35 = self.TensorProduct(ab35, uTable[4, :])r41 = self.TensorProduct(ab41, uTable[1, :])r42 = self.TensorProduct(ab42, uTable[2, :])r43 = self.TensorProduct(ab43, uTable[3, :])r44 = self.TensorProduct(ab44, uTable[4, :])r45 = self.TensorProduct(ab45, uTable[4, :])r51 = self.TensorProduct(ab51, uTable[2, :])r52 = self.TensorProduct(ab52, uTable[3, :])r53 = self.TensorProduct(ab53, uTable[4, :])r54 = self.TensorProduct(ab54, uTable[4, :])r55 = self.TensorProduct(ab55, uTable[4, :])rList = [r11, r12, r13, r14, r15,r21, r22, r23, r24, r25,r31, r32, r33, r34, r35,r41, r42, r43, r44, r45,r51, r52, r53, r54, r55]for rr in rList:R = self.MatricUnion(R, rr)if Err <= -40:e = np.array([1, 0, 0, 0, 0])if Err > -40 and Err <= -20:e = np.array([-1 / 20 * Err - 1, 1 / 20 * Err + 2, 0, 0, 0])if Err > -20 and Err <= 0:e = np.array([0, -1 / 20 * Err, 1 / 20 * Err + 1, 0, 0])if Err > 0 and Err <= 20:e = np.array([0, 0, -1 / 20 * Err + 1, 1 / 20 * Err, 0])if Err > 20 and Err <= 40:e = np.array([0, 0, 0, -1 / 20 * Err + 2, 1 / 20 * Err - 1])if Err > 40:e = np.array([0, 0, 0, 0, 1])if ErrK <= -40:ec = np.array([1, 0, 0, 0, 0])if ErrK > -40 and ErrK <= -20:ec = np.array([-1 / 20 * ErrK - 1, 1 / 20 * ErrK + 2, 0, 0, 0])if ErrK > -20 and ErrK <= 0:ec = np.array([0, -1 / 20 * ErrK, 1 / 20 * ErrK + 1, 0, 0])if ErrK > 0 and ErrK <= 20:ec = np.array([0, 0, -1 / 20 * ErrK + 1, 1 / 20 * ErrK, 0])if ErrK > 20 and ErrK <= 40:ec = np.array([0, 0, 0, -1 / 20 * ErrK + 2, 1 / 20 * ErrK - 1])if Err > 40:ec = np.array([0, 0, 0, 0, 1])u = self.Compose(self.ReshapeToRow(self.TensorProduct(e, ec)), R)u = u.squeeze()output = (0 * u[0] + 4 * u[1] + 8 * u[2] + 12 * u[3] + 16 * u[4]) / u.sum()output = np.clip(output, self.outputLast - 4, self.outputLast + 4)output = np.clip(output, self.outputMin, self.outputMax)self.outputLast = outputreturn outputdef CalDeltaErr(self, y):x = np.linspace(0, 9, num=10)n = x.shape[0]sum_xy = 0sum_x = 0sum_y = 0sum_xx = 0for i in range(n):sum_xy += x[i] * y[i]sum_x += x[i]sum_y += y[i]sum_xx += x[i] ** 2print(sum_xy, sum_x, sum_y, sum_xx)k = (n * sum_xy - sum_x * sum_y) / (n * sum_xx - sum_x ** 2)b = sum_y / n - k * sum_x / nreturn kif __name__ == '__main__':# 测试fc = FuzzyController(0)u = fc.Output(0, 0)print(u)print('done')
的控制域为、模糊集和隶属度函数按下图所示
的控制域、模糊集和隶属度函数与e相同。因此，e和
的模糊矩阵为
eTable = [[1, 0, 0, 0, 0],[0, 1, 0, 0, 0],[0, 0, 1, 0, 0],[0, 0, 0, 1, 0],[0, 0, 0, 0, 1]]
eTable = np.array(eTable)
ecTable = [[1, 0, 0, 0, 0],[0, 1, 0, 0, 0],[0, 0, 1, 0, 0],[0, 0, 0, 1, 0],[0, 0, 0, 0, 1]]
ecTable = np.array(ecTable)u的控制域、模糊集和隶属度函数如下图所示
u的模糊矩阵如下图所示
uTable = [[1, 0, 0, 0, 0],[0, 1, 0, 0, 0],[0, 0, 1, 0, 0],[0, 0, 0, 1, 0],[0, 0, 0, 0, 1]]
uTable = np.array(uTable)定义控制规则矩阵
代码实现如下
## 控制规则
a1 = eTable[0, :]
a2 = eTable[1, :]
a3 = eTable[2, :]
a4 = eTable[3, :]
a5 = eTable[4, :]b1 = ecTable[0, :]
b2 = ecTable[1, :]
b3 = ecTable[2, :]
b4 = ecTable[3, :]
b5 = ecTable[4, :]ab11 = self.ReshapeToCol(self.TensorProduct(a1, b1))
ab12 = self.ReshapeToCol(self.TensorProduct(a1, b2))
ab13 = self.ReshapeToCol(self.TensorProduct(a1, b3))
ab14 = self.ReshapeToCol(self.TensorProduct(a1, b4))
ab15 = self.ReshapeToCol(self.TensorProduct(a1, b5))
ab21 = self.ReshapeToCol(self.TensorProduct(a2, b1))
ab22 = self.ReshapeToCol(self.TensorProduct(a2, b2))
ab23 = self.ReshapeToCol(self.TensorProduct(a2, b3))
ab24 = self.ReshapeToCol(self.TensorProduct(a2, b4))
ab25 = self.ReshapeToCol(self.TensorProduct(a2, b5))
ab31 = self.ReshapeToCol(self.TensorProduct(a3, b1))
ab32 = self.ReshapeToCol(self.TensorProduct(a3, b2))
ab33 = self.ReshapeToCol(self.TensorProduct(a3, b3))
ab34 = self.ReshapeToCol(self.TensorProduct(a3, b4))
ab35 = self.ReshapeToCol(self.TensorProduct(a3, b5))
ab41 = self.ReshapeToCol(self.TensorProduct(a4, b1))
ab42 = self.ReshapeToCol(self.TensorProduct(a4, b2))
ab43 = self.ReshapeToCol(self.TensorProduct(a4, b3))
ab44 = self.ReshapeToCol(self.TensorProduct(a4, b4))
ab45 = self.ReshapeToCol(self.TensorProduct(a4, b5))
ab51 = self.ReshapeToCol(self.TensorProduct(a5, b1))
ab52 = self.ReshapeToCol(self.TensorProduct(a5, b2))
ab53 = self.ReshapeToCol(self.TensorProduct(a5, b3))
ab54 = self.ReshapeToCol(self.TensorProduct(a5, b4))
ab55 = self.ReshapeToCol(self.TensorProduct(a5, b5))r11 = self.TensorProduct(ab11, uTable[0, :])
r12 = self.TensorProduct(ab12, uTable[0, :])
r13 = self.TensorProduct(ab13, uTable[0, :])
r14 = self.TensorProduct(ab14, uTable[1, :])
r15 = self.TensorProduct(ab15, uTable[3, :])
r21 = self.TensorProduct(ab21, uTable[0, :])
r22 = self.TensorProduct(ab22, uTable[0, :])
r23 = self.TensorProduct(ab23, uTable[1, :])
r24 = self.TensorProduct(ab24, uTable[2, :])
r25 = self.TensorProduct(ab25, uTable[3, :])
r31 = self.TensorProduct(ab31, uTable[0, :])
r32 = self.TensorProduct(ab32, uTable[1, :])
r33 = self.TensorProduct(ab33, uTable[2, :])
r34 = self.TensorProduct(ab34, uTable[3, :])
r35 = self.TensorProduct(ab35, uTable[4, :])
r41 = self.TensorProduct(ab41, uTable[1, :])
r42 = self.TensorProduct(ab42, uTable[2, :])
r43 = self.TensorProduct(ab43, uTable[3, :])
r44 = self.TensorProduct(ab44, uTable[4, :])
r45 = self.TensorProduct(ab45, uTable[4, :])
r51 = self.TensorProduct(ab51, uTable[2, :])
r52 = self.TensorProduct(ab52, uTable[3, :])
r53 = self.TensorProduct(ab53, uTable[4, :])
r54 = self.TensorProduct(ab54, uTable[4, :])
r55 = self.TensorProduct(ab55, uTable[4, :])rList = [r11, r12, r13, r14, r15,r21, r22, r23, r24, r25,r31, r32, r33, r34, r35,r41, r42, r43, r44, r45,r51, r52, r53, r54, r55]
for rr in rList:R = self.MatricUnion(R, rr)拿出e为NB时e控制矩阵的对应行a1和为NB时
控制矩阵b1，进行张量乘法（注意顺序），对应函数为
def TensorProduct(self, a, b):将张量乘法的结果转为列向量（25*1）
def ReshapeToCol(self, a):e为NB为NB时，u为NB，因此，将上述所得列向量与u控制矩阵的NB对应行做张量乘法得到25行5列的矩阵r11，将规则矩阵的25个元素对应的r11到r55求出，25个矩阵取并集
def MatricUnion(self, a, b):得到规则矩阵R（25*5）
假设当e为-10，为20时，计算各自的隶属度矩阵（两个1*5的列向量），将e、
隶属度矩阵做张量乘法（注意顺序与上面对应），再转为行向量（1*25），与R矩阵（25*5）做矩阵合成运算
def Compose(self, a, b):得到u的隶属度矩阵（1*5），进行反模糊化，得到输出量
编写主程序main.py 运行仿真，查看控制器效果
from EPmodel import EPmodel
from fuzzy_control import FuzzyController
import numpy as np
import matplotlib.pyplot as plt
import pickle
import timeepm = EPmodel()
pd = 200
v = 30
n = 0
p = 300
fc = FuzzyController(n)
pCache = np.ones(10) * p
pIndex = 0
balanced = 0
nList = []
vList = []
pList = []
for i in range(1000):# fuzzy controlErrK = fc.CalDeltaErr(pCache)Err = p - pdn = fc.Output(Err=Err, ErrK=ErrK)p_ = epm.Pressure(p, v/60, n/60)pCache[pIndex] = p_pIndex += 1pIndex = pIndex % len(pCache)p = pCache.mean()print(v, p)pList.append(p_)nList.append(n)plt.figure('1')
plt.plot(pList)
plt.figure('2')
plt.plot(nList)plt.show()
本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！