使用python进行模糊控制模拟仿真

最新推荐文章于 2024-08-11 03:17:34 发布

河北一帆

最新推荐文章于 2024-08-11 03:17:34 发布

阅读量3k

点赞数 4

文章标签：模糊控制器 EPmodel 压力模拟张量乘法模糊逻辑

本文链接：https://blog.csdn.net/qq_41816368/article/details/109705906

版权

输入量为和 $\Delta e$ ，单输出

构建一个模拟环境EPmodel.py，一个输入一个输出

import numpy as np
import matplotlib.pyplot as plt
import math
import random


Pa = 101

class EPmodel:
    def __init__(self):
        self.a = 973
        self.n = 0.955
        self.eta = 0.775
        self.diameter = 6.280  
        self.A = 1/4 * math.pi * self.diameter * self.diameter
        self.L = 900  
        self.Vc = self.A * self.L * 0.599
        self.As = 0.953  
        self.h = 960  

    def Pressure(self, p, v, n):
        K1 = p / Pa
        K2 = self.a * math.pow(K1, self.n) / self.Vc

        K3 = self.A * v
        K4 = (self.As * self.h * self.eta * n) / (2 * math.pi)
        K5 = K3 - K4
        deltaP = K2 * K5

        p_ = p + deltaP

        p__ = np.random.normal(loc=p_, scale=2.0, size=1)

        return p__[0]


if __name__ == '__main__':
    #测试模拟环境
    epm = EPmodel()
    p = 120
    v = 40
    n = 15
    v = v / 60
    n = n / 60

    pList = []
    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 np


class FuzzyController:
    def __init__(self, n0):
        self.outputLast = n0
        self.outputMin = 0
        self.outputMax = 16

    def 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 x
        for i in range(rows):
            for j in range(cols):
                c[i, j] = ChooseSmall(a[i], b[j])

        return c

    def ReshapeToCol(self, a):
        rows = a.shape[0]
        cols = a.shape[1]
        length = rows * cols
        b = np.ones(length)
        for i in range(rows):
            r = i*rows
            c = i*rows+cols
            b[r:c] = a[i, :]

        return b[..., np.newaxis]

    def ReshapeToRow(self, a):
        rows = a.shape[0]
        cols = a.shape[1]
        length = rows * cols
        b = np.ones(length)
        for i in range(rows):
            r = i * rows
            c = i * rows + cols
            b[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 c

    def VectorIntersection(self, a, b):
        ##两个向量取交集（取小数）
        ChooseSmall = lambda x, y: y if x > y else x
        c = np.ones(a.shape[0])
        for i in range(c.shape[0]):
            c[i]= ChooseSmall(a[i], b[i])

        return c

    def 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 c

    def Output(self, Err, ErrK):
        ## Err是当前值-目标值 x - xd
        # -40, -20, 0, 20, 40
        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)
        # -40, -20, 0, 20, 40
        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)
        # -4, -2, 0, 2, 4
        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)
        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 = output

        return output

    def CalDeltaErr(self, y):
        x = np.linspace(0, 9, num=10)
        n = x.shape[0]
        sum_xy = 0
        sum_x = 0
        sum_y = 0
        sum_xx = 0
        for i in range(n):
            sum_xy += x[i] * y[i]
            sum_x += x[i]
            sum_y += y[i]
            sum_xx += x[i] ** 2

        print(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 / n

        return k


if __name__ == '__main__':
    # 测试
    fc = FuzzyController(0)
    u = fc.Output(0, 0)
    print(u)

    print('done')

的控制域为、模糊集和隶属度函数按下图所示

$\Delta e$ 的控制域、模糊集和隶属度函数与e相同。因此，e和 $\Delta 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和 $\Delta e$ 为NB时 $\Delta e$ 控制矩阵b1，进行张量乘法（注意顺序），对应函数为

def TensorProduct(self, a, b):

将张量乘法的结果转为列向量（25*1）

def ReshapeToCol(self, a):

e为NB $\Delta e$ 为NB时，u为NB，因此，将上述所得列向量与u控制矩阵的NB对应行做张量乘法得到25行5列的矩阵r11，将规则矩阵的25个元素对应的r11到r55求出，25个矩阵取并集

def MatricUnion(self, a, b):

得到规则矩阵R（25*5）

假设当e为-10， $\Delta e$ 为20时，计算各自的隶属度矩阵（两个1*5的列向量），将e、 $\Delta 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 time


epm = 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 control
    ErrK = fc.CalDeltaErr(pCache)
    Err = p - pd
    n = fc.Output(Err=Err, ErrK=ErrK)

    p_ = epm.Pressure(p, v/60, n/60)
    pCache[pIndex] = p_
    pIndex += 1
    pIndex = 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()