美文网首页Python
Python实现《Robust output feedback

Python实现《Robust output feedback

作者: Python解决方案 | 来源:发表于2019-12-27 09:20 被阅读0次

    Robust output feedback fault-tolerant control of non-linear multi-agent systems based on wavelet neural networks摘要如下:
    A robust output feedback active fault-tolerant leader-following controller for a class of non-linear multi-agent systems is presented. It is assumed that the states of the followers are not available; therefore, a local observer is constructed to estimate the states of each agent. In addition, the non-linear dynamics of agents may include uncertainties and the control input
    of the leader dynamics is unknown to all followers. Moreover, taking advantage of wavelet neural networks (WNNs), an online fault estimation scheme is developed which can effectively approximate the unknown actuator faults. The proposed decentralised observer-based robust cooperative controller is capable of compensating for the effects of unknown time-varying additive actuator faults, the model uncertainties, and the unknown input of the leader simultaneously. The stability analysis and convergence results that guarantee boundedness of all closed-loop signals are investigated via Lyapunov's direct method. To demonstrate the effectiveness of the proposed approach, a network of single-link manipulators is studied. As the results verify, the proposed WNN-based fault estimation scheme can properly approximate the unknown actuator faults, which results in efficient compensation in the fault-tolerant control design to achieve cooperative tracking objectives.
    仿真实验代码如下:

    #导入需要用到的库
    #-------------------------------------------------------------------------------------------------------------------------
    import math
    import numpy as np
    import matplotlib.pyplot as plt
    #-------------------------------------------------------------------------------------------------------------------------
    
    #定义前向控制系统相关参数
    #-------------------------------------------------------------------------------------------------------------------------
    t=0
    c = 4
    rho = 0.1
    epsilon = 0.01
    A = np.mat([[0,1,0,0],[-48.6,-1.25,48.6,0],[0,0,0,10],[1.95,0,-1.95,0]])#(4x4)
    B = np.mat([0,21.6,0,0]).T#(4x1)
    C = np.mat([[1,0,0,0],[0,1,0,0]])#(2x4)
    D = np.mat([0,1,0,0]).T#(4x1)
    E = np.mat([0,12.5,0,0]).T#(4x1)
    K = np.mat([-2.2297,-0.4803,1.3470,-2.3819])#(1x4)
    An =np.mat([[0,1,0,0,1],[1,0,1,0,0],[0,1,0,1,0],[0,0,1,0,1],[1,0,0,1,0]])#(5x5)#跟随者连通矩阵
    Bn = np.mat([[1,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0]])#(5x5)#领导者连通矩阵
    F1 = np.mat([[1.8275,0.9931],[0.9931,13.7082],[1.8124,2.8758],[0.1553,0.5663]])#(4x2)
    F2 = np.mat([3.9416,1.0812])
    P1 = np.mat([[4.5785,0.3153,-2.5338,1.2168],
                [0.3153,0.0865,0.0,0.0],
                [-2.5338,0.0,2.4053,-1.4706],
                [1.2168,0.0,-1.4706,4.6813]])
    #-------------------------------------------------------------------------------------------------------------------------
    
    #各类状态及输出初始化
    #-------------------------------------------------------------------------------------------------------------------------
    X_leader=np.mat(np.random.rand(4,1))#领导者状态初始化,可进行随机初始化
    Y_leader=np.mat(np.random.rand(2,1))#领导者输出零初始化
    
    X=np.mat(np.random.rand(4,5))#跟随者状态随机初始化,第一列表示第一个跟随者的状态
    Y=np.mat(np.random.rand(2,5))#跟随者输出随机初始化,第一列表示第一个跟随者的输出
    X_hat=np.mat(np.zeros((4,5)))#估计器状态零初始化,第一列表示第一个跟随者的估计状态
    Y_hat=np.mat(np.zeros((2,5)))#估计器输出零初始化,第一列表示第一个跟随者的估计输出
    X_wav=np.mat(np.zeros((4,5)))#状态偏差零初始化,第一列表示第一个跟随者的状态状态
    Y_wav=np.mat(np.zeros((2,5)))#输出偏差零初始化,第一列表示第一个跟随者的输出偏差
    #-------------------------------------------------------------------------------------------------------------------------
    
    #定义小波神经网络的相关参数
    #-------------------------------------------------------------------------------------------------------------------------
    z=5#小波网络的个数
    m=2#输出的分量个数m
    p=1#网络输出层个数1
    eta1=2#超参数eta1
    eta2=1#超参数eta2
    eta3=1#超参数eta3
    delta1=1#公式(14)求积部分
    delta2=0#公式(14)求和部分
    phi_mu=0#偏导数的中间量phi_mu
    phi_sig=0#偏导数的中间量phi_sig
    
    T_hat=np.mat([0,0,0,0,0])#传动误差估计值零初始化
    G_mu = np.zeros((5,15,5))#见公式(20)
    G_sig=np.zeros((5,15,5))#见公式(20)
    W = np.zeros((5,5))#小波网络权重参数初始化
    M = np.mat(np.zeros((3,5)))#小波网络输入初始化
    W_hat = np.mat(np.zeros((5,5)))#权重估计参数(5x1x5),W_i(5x1)
    mu=np.mat(np.zeros((15,5)))#膨胀参数(zx(m+p)x5)
    mu_hat=np.mat(np.zeros((15,5)))#膨胀参数(zx(m+p)x5)
    sig=np.mat(np.zeros((15,5)))#转换参数(zx(m+p)x5)
    sig_hat=np.mat(np.zeros((15,5)))#转换参数(zx(m+p)x5)
    Phi = np.mat(np.zeros((5,5)))#小波方程初始化Phi_i(5x1)
    Phi_hat=np.mat(np.zeros((5,5)))#小波方程估计初始化
    #-------------------------------------------------------------------------------------------------------------------------
    
    #定义保存各种中间量的列表
    #-------------------------------------------------------------------------------------------------------------------------
    #保存领导者的状态和输出
    X_leader_list=[]
    Y_leader_list=[]
    #保存跟随者的状态和输出
    X_list=[]
    Y_list=[]
    #保存跟随着的估计状态和估计输出
    X_hat_list=[]
    Y_hat_list=[]
    #保存跟随者估计状态偏差和估计输出偏差
    X_wav_list=[]
    Y_wav_list=[]
    #保存每个跟随者的偏导数的每一个分量,共15个
    phi_mu_list=[]
    phi_sig_list=[]
    #保存传动错误和估计传动错误
    T_list=[]
    T_hat_list=[]
    #保存跟随着者与领导者的状态偏差以及控制器的中间量
    E_list=[]
    U_list=[]
    #-------------------------------------------------------------------------------------------------------------------------
    
    #定义T_bag函数,只有第二个跟随者和第四个跟随者存在actuator fault,并且只出现在特定的时间段(1x5)
    #-------------------------------------------------------------------------------------------------------------------------
    def T_bag(t):
        if (t>=0 and t<5) or t>=30:
            T=np.mat([0,0,0,0,0])
        elif t>=5 and t<8:
            T=np.mat([0,0.4*math.sin(t)*math.cos(0.3*t),0,0,0])
        else:
            T=np.mat([0,0.4*math.sin(t)*math.cos(0.3*t),0,0.3*math.sin(0.4*t),0])
        return T
    #-------------------------------------------------------------------------------------------------------------------------
    
    #学习率lr
    lr=0.002
    
    #在时长为40的范围内进行更新迭代
    #-------------------------------------------------------------------------------------------------------------------------
    for t in np.arange(0,40,lr):
    
    #0.领导者的状态更新
    #-------------------------------------------------------------------------------------------------------------------------
        part1=np.dot(A,X_leader)#A*x1   (4x1)
        part2=np.mat([0,0,0,-0.333*math.sin(X_leader[2,0])]).T#g(x1)   (4x1)
        part3=np.dot(B,0.1*math.sin(0.1*t))#B*u1  (4x1)
        dx_leader=part1+part2+part3
        X_leader=X_leader+dx_leader*lr
        Y_leader=np.dot(C,X_leader)
        
        X_leader_list.append(X_leader)
        Y_leader_list.append(Y_leader)
    #-------------------------------------------------------------------------------------------------------------------------
    
    #1.计算U,U=Un+Ur+Uac
    #-------------------------------------------------------------------------------------------------------------------------        
        Un=[]
        Ur=[]
        Uac=[]
        U=[]#跟随者控制器初始化
        for i in range(5):
            inner_sum=np.mat([0.,0.,0.,0.]).T
            for j in range(5):
                inner_sum+=An[j,i]*(X_hat[:,i]-X_hat[:,j])#(4x1)
            Un.append(c*np.dot(K,(inner_sum+Bn[i,i]*(X_hat[:,i]-X_leader))))#(1x1 scalar)
            Ur.append(rho*np.tanh(np.dot(K,(inner_sum+Bn[i,i]*(X_hat[:,i]-X_leader)))/epsilon))#(1x1 scalar)
            Uac.append(-np.dot(B.I,E*T_hat[0,i]))#(1x1 scalar)
        U=np.mat(np.array(Un))+np.mat(np.array(Ur))+np.mat(np.array(Uac))#(1x5)m每一列代表一个跟随者的控制器
        #将U存进U_list
        U_list.append(U)
    #-------------------------------------------------------------------------------------------------------------------------
    
    #2.计算dX dX=AX+g(x)+delta_f(x)+Bu+ET+Dw
    #-------------------------------------------------------------------------------------------------------------------------
        p1=np.dot(A,X)#(4x5)
        mat0=np.mat([0,0,0,0,0])
        p2=np.concatenate((mat0,mat0,mat0,-0.333*np.sin(X[2,:])),axis=0)#(4x5)
        p3=np.concatenate((mat0,0.1*np.cos(X[1,:]),mat0,mat0),axis=0)#(4x5)
        p4=np.dot(B,U)#(4x5)
        p5=np.dot(E,T_bag(t))#(4x5)
        p6=np.dot(D,0.1*math.sin(10*t))#(4x1)
        dX=(p1+p2+p3+p4+p5+p6)#(4X5)
    #-------------------------------------------------------------------------------------------------------------------------
    
    #3.计算跟随者每个时刻更新后的状态和输出
    #-------------------------------------------------------------------------------------------------------------------------
        X=X+dX*lr#(4x5)
        Y=np.dot(C,X)#(2x5)
        X_list.append(X)#(4000x4x5)
        Y_list.append(Y)#(4000x2x5)
    #-------------------------------------------------------------------------------------------------------------------------
        
    #4.计算dX_hat
    #-------------------------------------------------------------------------------------------------------------------------
        p_hat1=np.dot(A,X_hat)#(4x5)
        p_hat2=np.concatenate((mat0,mat0,mat0,-0.333*np.sin(X_hat[2,:])),axis=0)#(4x5)
        p_hat3=np.dot(B,U)#(4x5)
        p_hat4=np.dot(E,T_hat)#(4X5)
        p_hat5=np.dot(F1,(Y-Y_hat))#(4x2)*(2*5)-->(4x5)
        ppp=np.dot(C.T,np.sign(0.015*(Y-Y_hat)))
        p_hat6=0.015*np.dot(P1.I,ppp)#(4x4)*(4x2)*(2x5)-->(4x5)
        dX_hat=p_hat1+p_hat2+p_hat3+p_hat4+p_hat5+p_hat6
    #-------------------------------------------------------------------------------------------------------------------------
        
    #5.计算观测器每个时刻更新后的状态和输出并保存,计算估计器和执行器的状态偏差和输出偏差并保存
    #-------------------------------------------------------------------------------------------------------------------------
        X_hat=X_hat+dX_hat*lr#(4x5)
        Y_hat=np.dot(C,X_hat)#(2x5)
        X_wav=X-X_hat#(4x5)
        Y_wav=Y-Y_hat#(2x5)
        
        X_hat_list.append(X_hat)#(40000x4x5)
        Y_hat_list.append(Y_hat)#(40000x2x5)
        X_wav_list.append(X_wav)#(40000x4x5)
        Y_wav_list.append(Y_wav)#(40000x2x5)
    #-------------------------------------------------------------------------------------------------------------------------
    
    #6.求解FTC中关键组成部分T_hat
    #-------------------------------------------------------------------------------------------------------------------------
        #(1).M由两部分组成,见公式(12)下文字
        M=np.vstack((Y_wav,T_hat))#(3x5)#小波函数输入
        
        #(2).求解Phi、W_hat、mu_hat(该程序用mu代替)、sig_hat(该程序用sig代替)
        for i in range(5):
            for j in range(5):
                delta1=1
                delta2=0
                for k in range(3):
                    delta1=delta1*(1-mu[j*3+k,i]*mu[j*3+k,i]*M[k,i]*M[k,i])#公式(14)求积部分
                    delta2=delta2+mu[j*3+k,i]*mu[j*3+k,i]*(M[k,i]-sig[j*3+k,i])*(M[k,i]-sig[j*3+k,i])#公式(14)求和部分
                Phi[j,i]=delta1*np.exp(-0.5*delta2)#公式(14)整体
                for k in range(3):
                    #公式(21)里的Omiga_phi_ij/Omiga_mu_ijk
                    phi_mu=((-2*M[k,i]*M[k,i])/(1-mu[j*3+k,i]*mu[j*3+k,i]*M[k,i]*M[k,i])
                            -np.power((M[k,i]-sig[j*3+k,i]),2))*mu[j*3+k,i]*delta1*np.exp(-0.5*delta2)
                    #公式(21)里的Omiga_phi_ij/Omiga_sig_ijk
                    phi_sig=delta1*np.exp(-0.5*delta2)*mu[j*3+k,i]*mu[j*3+k,i]*(M[k,i]-sig[j*3+k,i])
                    #先分别保存在各自list中
                    phi_mu_list.append(phi_mu)
                    phi_sig_list.append(phi_sig)
                    
            #将list中的值分解到各个矩阵中,先组成公式(21),再组成(20)
            for j in range(5):
                G_mu[i][j*3:(j+1)*3][:,j]=phi_mu_list[j*3:(j+1)*3]
                G_sig[i][j*3:(j+1)*3][:,j]=phi_sig_list[j*3:(j+1)*3]
            
            #通过公式(25)更新W_hat、mu、sig
            dWi_hat=eta1*(Phi[:,i]+np.dot(G_mu[i].T,mu[:,i])-np.dot(G_sig[i].T,sig[:,i]))*((np.dot(F2,Y_wav[:,i])).T)
            W_hat[:,i]=W_hat[:,i]+dWi_hat*lr
            mu[:,i]=eta2*np.dot(G_mu[i],W_hat[:,i]*np.dot(F2,Y_wav[:,i]))
            sig[:,i]=eta3*np.dot(G_sig[i],W_hat[:,i]*np.dot(F2,Y_wav[:,i]))
            
        #(3).得到M、Phi后,即可更新T_hat
        T_hat=np.mat(np.diag(np.dot(W_hat.T,Phi)))
        T_hat_list.append(T_hat)
        
        #(4).保存每个真正的T
        T_list.append(T_bag(t))
    #-------------------------------------------------------------------------------------------------------------------------
    

    画图:

    num=int(30/lr)
    t=np.linspace(0,30,num=num)
    fig = plt.figure(figsize=(16,4))
    
    ax1, ax2 = fig.subplots(1, 2, sharey=True)
    ax1.plot(t,np.array(T_list)[:num,:2,1],'b-')
    ax1.plot(t,np.array(T_hat_list)[:num,:2,1],'r--')
    ax1.legend(('T2', 'Estimated T2'), loc='upper right') 
    ax1.set_xlabel('Time(s)')
    ax1.set_title('Fault Estimation T2')
    ax1.grid()
    
    ax2.plot(t,np.array(T_list)[:num,:2,3],'b-')
    ax2.plot(t,np.array(T_hat_list)[:num,:2,3],'r--')
    plt.legend(('T4', 'Estimated T4'), loc='upper right')
    ax2.set_xlabel('Time(s)')
    ax2.set_title('Fault Estimation T4')
    ax2.grid()
    plt.show()
    
    Fault Estimation
    num=int(40/lr)
    t=np.linspace(0,40,num=40/lr)
    fig = plt.figure(figsize=(16,8))
    axes = fig.subplots(2, 2, sharey=True)
    
    axes[0,0].plot(t,np.array(X_wav_list)[:num,0,:5])
    axes[0,0].legend(("Agent1", "Agent2","Agent3","Agent4","Agent5"), loc='upper right') 
    axes[0,0].set_ylabel('1st State')
    axes[0,0].set_title('State Estimation Errors')
    axes[0,0].grid()
    
    axes[0,1].plot(t,np.array(X_wav_list)[:num,1,:5])
    axes[0,1].legend(("Agent1", "Agent2","Agent3","Agent4","Agent5"), loc='upper right') 
    axes[0,1].set_ylabel('2nd State')
    axes[0,1].set_title('State Estimation Errors')
    axes[0,1].grid()
    
    axes[1,0].plot(t,np.array(X_wav_list)[:num,2,:5])
    axes[1,0].legend(("Agent1", "Agent2","Agent3","Agent4","Agent5"), loc='upper right') 
    axes[1,0].set_xlabel('Time(s)')
    axes[1,0].set_ylabel('3rd State')
    axes[1,0].grid()
    
    axes[1,1].plot(t,np.array(X_wav_list)[:num,3,:5])
    axes[1,1].legend(("Agent1", "Agent2","Agent3","Agent4","Agent5"), loc='upper right') 
    axes[1,1].set_xlabel('Time(s)')
    axes[1,1].set_ylabel('4th State')
    axes[1,1].grid()
    plt.show()
    
    State Estimation Errors
    num=int(40/lr)
    t=np.linspace(0,40,num=40/lr)
    fig = plt.figure(figsize=(16,8))
    axes = fig.subplots(2, 2, sharey=True)
    
    axes[0,0].plot(t,(np.array(X_list)-np.array(X_leader_list))[:num,0,:5])
    axes[0,0].legend(("Agent1", "Agent2","Agent3","Agent4","Agent5"), loc='upper right') 
    axes[0,0].set_ylabel('1st State')
    axes[0,0].set_title('Leader-Following Tracking Errors')
    axes[0,0].grid()
    
    axes[0,1].plot(t,(np.array(X_list)-np.array(X_leader_list))[:num,1,:5])
    axes[0,1].legend(("Agent1", "Agent2","Agent3","Agent4","Agent5"), loc='upper right')
    axes[0,1].set_ylabel('2nd State')
    axes[0,1].set_title('Leader-Following Tracking Errors')
    axes[0,1].grid()
    
    axes[1,0].plot(t,(np.array(X_list)-np.array(X_leader_list))[:num,2,:5])
    axes[1,0].legend(("Agent1", "Agent2","Agent3","Agent4","Agent5"), loc='upper right') 
    axes[1,0].set_xlabel('Time(s)')
    axes[1,0].set_ylabel('3rd State')
    axes[1,0].grid()
    
    axes[1,1].plot(t,(np.array(X_list)-np.array(X_leader_list))[:num,3,:5])
    axes[1,1].legend(("Agent1", "Agent2","Agent3","Agent4","Agent5"), loc='upper right') 
    axes[1,1].set_xlabel('Time(s)')
    axes[1,1].set_ylabel('4th State')
    axes[1,1].grid()
    plt.show()
    
    Leader-Following Tracking Errors
    num=int(40/lr)
    t=np.linspace(0,40,num=40/lr)
    plt.figure(figsize=(8,6))
    plt.plot(t,np.array(U_list)[:num,0,:5])
    plt.legend(("Agent1", "Agent2","Agent3","Agent4","Agent5"), loc='upper right') 
    plt.xlabel('Time(s)')
    plt.title('Control Signals')  
    plt.grid()
    plt.show()
    
    Control Signals

    相关文章

      网友评论

        本文标题:Python实现《Robust output feedback

        本文链接:https://www.haomeiwen.com/subject/roproctx.html