BP网络的代码：

function [sys,x0,str,ts]=BPwangluo(t,x,u,flag)

switch flag

case 0

[sys,x0,str,ts]=mdlInitializeSizes;

case 2

sys=mdlUpdates(x,u);

case 3

sys=mdlOutputs(t,x,u);

case {1,4,9}

sys=[];

otherwise

error(['unhandled flag=',num2str(flag)]);%异常处理

end

function[sys,x0,str,ts]=mdlInitializeSizes

sizes=simsizes;%用于设置模块参数的结构体用simsizes来生成

sizes.NumContStates=0;%模块连续状态变量的个数

sizes.NumDiscStates=6;%模块离散状态变量的个数

sizes.NumOutputs=3;%模块输出变量的个数

sizes.NumInputs=14;%模块输入变量的个数

sizes.DirFeedthrough=1;%模块是否存在直接贯通，1表示存在直接贯通，若为0，则mdlOutputs函数里不能有u

sizes.NumSampleTimes=1;%模块的采样时间个数,至少是一个

sys=simsizes(sizes);%设置完后赋给sys输出

x0=zeros(6,1);%系统状态变量设置

str=[];

ts=[1 0];%采样周期设为0表示是连续系统，

function sys=mdlUpdates(x,u)

x1=[u(5);x(2)+u(5);(u(5)-u(6))];%3个状态量(电机的偏差、偏差和以及偏差变化量)，u(5)是偏差，u(4)是上一次的偏差，x(2)是之前的偏差和

x2=[u(7);x(5)+u(7);(u(7)-u(8))];%3个状态量(电机的偏差、偏差和以及偏差变化量)，u(7)是偏差，u(8)是上一次的偏差，x(5)则是之前的偏差和

sys=[x1(1);x1(2);x1(3);x2(1);x2(2);x2(3)];

function sys=mdlOutputs(t,x,u)

xite=0.2;

alfa=0.05;

IN=6;H=8;OUT=3;

wi=rand(8,6);%产生一个的随机8*6数矩阵，随机数在(0，1)区间

wi_1=wi;wi_2=wi;wi_3=wi;wi_4=wi;wi_5=wi;wi_6=wi;

wo=rand(3,8);%产生一个的随机3*8数矩阵，随机数在(0，1)区间

wo_1=wo;wo_2=wo;wo_3=wo;wo_4=wo;wo_5=wo;wo_6=wo;wo_7=wo;wo_8=wo;

Oh=zeros(1,8);%产生一个1*8的零矩阵(行矩阵)

I=Oh;

xi=[u(9),u(5),u(3),u(10),u(7),u(11)];%神经网络的6个输入，电机2的期望值、误差以及实际值，电机1的期望值、误差以及实际值

emfa1=[x1(1);x1(2);x1(3)];emfa2=[x2(1);x2(2);x2(3)]%6个状态变量(两个电机的偏差、偏差和、偏差变化量)(两个3*1矩阵，列向量)

I=xi*wi';%隐层的输入(1*8的矩阵)

for j=1:1:8

Oh(j)=(exp(I(j))-exp(-I(j)))/(exp(I(j))+exp(-I(j)));%隐层的输出值(1*8矩阵)行矩阵

end

K1=Oh*wo';%输出层的输入(1*3矩阵)

for i=1:1:3

K(i)=exp(K1(i))/(exp(K1(i))+exp(-K1(i)));%得到输出层的输出：惩罚因子、学习率(1*3矩阵，行向量)

end

%%以下是权值调整

%隐含层至输出层的权值调整

dyu1=sign((u(3)-u(4))/(u(1)-u(2)+0.0001));

dyu2=sign((u(11)-u(12))/(u(13)-u(14)+0.0001))

for j=1:1:3

dK(j)=2/(exp(K1(j))+exp(-K1(j)))^2; %输出层的输出的一阶导

end

for i=1:1:3

delta31(i)=u(5)*dyu1*emfa1(i)*dK(i);  %输出层的delta1

delta32(i)=u(7)*dyu2*emfa2(i)*dK(i);  %输出层的delta2

delta3(i)= (delta31(i)+delta32(i))/2;  %输出层的delta的均值,3*1矩阵，列向量

end

for j=1:1:3

for i=1:1:8

d_wo=xite*delta3(j)*Oh(i)+alfa*(wo_1-wo_2);

end

end

wo=wo_1+d_wo;   %3*8矩阵

%以下是输入层至隐含层的权值调整

for i=1:1:8

dO(i)=4/(exp(I(i))+exp(-I(i)))^2;%(1*8矩阵)

end

segma=delta3'*wo; %(1*8矩阵，行向量)

delta2 = dO.*segma; %1*8矩阵

d_wi = xite*delta2'*xi+alfa*(wi_1-wi_2);

wi=wi_1+d_wi;

wo_6=wo_5;

wo_5=wo_4;

wo_4=wo_3;

wo_3=wo_2;

wo_2=wo_1;

wo_1=wo;   %储存输出层本次调整后的权值

wi_6=wi_5;

wi_5=wi_4;

wi_4=wi_3;

wi_3=wi_2;

wi_2=wi_1;

wi_1=wi;   %储存隐层本次调整后的权值

sys=[K(1),K(2),K(3)];