一、代码获取方式获取代码方式1：
完整的代码上传了我的资源。 【空中程序】基于matlab RBF神经网络控制卫星轨道和姿态的【含Matlab源代码的377期】
获取代码方式2：
订阅紫极神光博客3358www.Sina.com/，可在支付证书付费专栏上获得此代码。
评论：
订阅紫极神光博客3358www.Sina.com/，即可免费获得私信博主部的代码(付费专栏自订阅之日起3天有效)。
二、RBF个人资料1
1985年，Powell提出了多元插值的径向基函数(RBF )方法。 径向基函数是取只依赖于距原点的距离的值的实数函数。 也就是说，到(x (=x )或任意一点c的距离也可以。 c点称为中心点。 即，(x，c(=x-c )。 满足(x )=(x ) ) )特性的函数均称为径向基函数，标准上一般使用rxdqt距离(也称为欧式径向基函数) )，但也可以是其他距离函数。 最常用的径向基函数是tsdyc核函数，形式为k(|x-xc|| )=exp(|x-xc||2/(2*)2) )，x_c为核函数的中心，为函数的宽度参数，函数的
有效期
RBF神将网络是三层神经网络，包括输入层、隐藏层和输出层。 从输入空间到隐层空间的转换是非线性的，从隐层空间到输出层空间的转换是线性的。 流程图如下。
RBF网络的基本思想是以RBF为隐藏单元的“基础”构成隐藏层空间，允许输入向量直接映射到隐藏空间而无需加权连接。 RBF的中心点确定后，该映射关系也确定。 隐式层空间到输出空间的映射是线性的，即，网络的输出是隐式单元的输出的线性加权和，并且这里的权重是网络可调参数。 其中，隐式层的作用通过将向量映射到低维p到qpdmb的h，使得能够从低维的不可线性分离情形到qpdmb线性分离。 主要是核函数的思想。 这样，网络从输入到输出的映射是非线性的，并且网络的输出相对于可调参数是线性的。 网络权重可以直接由线性方程求解，大大提高学习速度，避免局部极小问题。
径向基神经网络的激活函数可以表示为。
其中，xp是第p个输入样本，ci是第I个中心点，h是隐含层中的节点数，且n是输出的样本数或分类数。 从放射状神经网络的结构中得到的网络的输出如下
当然，使用最凶暴的耳机乘坐的损失函数表示如下。
1什么是径向基函数
求解的参数有三个：基函数的中心、方差和从隐式层到输出层的权重。
(1)自组织中心学习方法：
第一步：不监控学习过程，求隐式层基函数的中心和方差
第二步：监督学习过程，求解隐含层到输出层之间的权重
首先，选择h个中心作为k均值聚类，用公式求出tsdyc核函数径向基的方差：
cmax是选定中心点之间的最大距离。
从抑制层到输出层之间的神经元连接权重可以用最凶恶的耳机乘法直接计算，即对于损耗函数，将w的偏导数求解为等于0，可以得到以下简并计算公式。
)2)直接计算法
抑制层神经元的中心选自随机输入样本，中心固定。 中心固定后，隐式层神经元的输出是已知的，这种神经网络的连接权可以通过求解线性方程组来确定。 应用于样本数据的分布明显具有代表性。
)3)有监控学习算法
训练样本集获得满足监测要求的网络中心和其他权重参数，经过误差修正学习过程，与BP网络学习原理一样采用梯度下降法。 因此，RBF也可以看作是BP神经网络的一种。
http://www.Sina.com/http://www.Sina.com /
BP神经网络的隐藏节点采用输入模式与权向量内积作为激活函数的自变量，激活函数采用Sigmoid函数。 由于各参数对BP网络的输出有同等的影响，BP神经网络是对非线性映射的全局逼近。
RBF神经网络的隐藏节点使用输入模式与中心向量之间的距离(例如欧式距离)作为函数的参数，使用径向基函数(例如Gaussian函数)作为激活函数。 神经元的输入越远离径向基函数的中心，神经元的激活程度越低(tsdyc函数)。 RBF神经网络的输出与部分参数相关，例如，一个wij值只影响一个yi的输出，因此RBF神经网络具有“局部映射”特性。
局部近似是指目标函数的近似只基于查询点附近的数据。 实际上，在径向基网络中，通常使用tsdyc径向基函数。 函数图像是双侧衰减的、径向对称的，并且只有当所选择的中心接近查询点(即输入数据)时，输入才具有真实的映射作用。 中心与查询点较远、欧式距离过大时，输出结果接近0，因此真实作用点是接近查询点的点，局部近似； 对BP网络目标函数的近似不仅涉及查询点附近的数据，还涉及所有数据。
ng>4.2 中间层数的区别BP神经网络可以有多个隐含层，但是RBF只有一个隐含层。
4.3 训练速度的区别
使用RBF的训练速度快，一方面是因为隐含层较少，另一方面，局部逼近可以简化计算量。对于一个输入x，只有部分神经元会有响应，其他的都近似为0，对应的w就不用调参了。
4.4 RBF网络是连续函数的最佳逼近，而BP网络不是。
三、部分源代码 clc;clear all;close all;mu = 3.986004418e14;al = 7400000; %参考星轨道半长轴(m)n = sqrt(mu/(al^3));%参考星平均角速度N = 2;T = 2*pi/n;t0 = 0;tf = t0 + N*T;x0=[0.3;-0.1;0.2;-0.2;0.12;0.17;-0.09];% xt0=[0.3;-0.1;0.2;-0.2;0.12;0.17;-0.09;-1.3;-3.05;1.22];% [t,x]=ode45('zitai_function',[t0:17:tf],x0);% [t,xt]=ode45('zitai_ctrl',[t0:17:tf],xt0); % [t,xq]=ode45('oula',[t0:17:tf],[-1.3;-3.05;1.22]);[t,x]=ode45('rbf_zitai',[t0:17:5000],x0);% c0=pi/T;% wrd=[-c0*sin(2*c0*t) 1.6*c0*sin(4*c0*t) 0.8*c0*sin(2*c0*t)];% wr=[xt(:,5)-wrd(1,1) xt(:,6)-wrd(2,1) xt(:,7)-wrd(3,1)];% fai=xq(:,1);sita=xq(:,2);pesi=xq(:,3);% for i=1:1118% q(i,1)=cos(pesi(i,:)/2)*sin(fai(i,:)/2)*cos(sita(i,:)/2)-sin(pesi(i,:)/2)*cos(fai(i,:)/2)*sin(sita(i,:)/2);% q(i,2)=cos(pesi(i,:)/2)*cos(fai(i,:)/2)*sin(sita(i,:)/2)+sin(pesi(i,:)/2)*sin(fai(i,:)/2)*cos(sita(i,:)/2);% q(i,3)=sin(pesi(i,:)/2)*cos(fai(i,:)/2)*cos(sita(i,:)/2)+cos(pesi(i,:)/2)*sin(fai(i,:)/2)*sin(sita(i,:)/2);% q(i,4)=cos(pesi(i,:)/2)*cos(fai(i,:)/2)*cos(sita(i,:)/2)-sin(pesi(i,:)/2)*sin(fai(i,:)/2)*sin(sita(i,:)/2);% endunction dx=rbf_micishedong(t,x)mu=398600.44e+9;rol=7400000;Re = 6378000;%地球赤道半径el = 0; %偏心率il = 30*pi/180; %轨道倾角Ml = 90*pi/180; %平近点角nl = sqrt(mu/rol^3);wl = 0*pi/180; %近地点辐角raanl = 100*pi/180; %升交点赤经fl = 160.0513*pi/180;%真近点角ul = wl+fl;%真纬度角J2 = 1.08263e-3;A=[cos(ul)*cos(raanl)-sin(ul)*cos(il)*sin(raanl) cos(ul)*sin(raanl)+sin(ul)*cos(il)*cos(raanl) sin(ul)*sin(il); -sin(ul)*cos(raanl)-cos(ul)*cos(il)*sin(raanl) -sin(ul)*sin(raanl)+cos(ul)*cos(il)*cos(raanl) cos(ul)*sin(il); sin(il)*sin(raanl) -sin(il)*cos(raanl) cos(il)];%地心惯性坐标系到参考航天器的坐标转换矩阵fal1=sqrt(mu*rol*(1-el^2))*(1+2*el*cos(Ml)+el^2*0.5*(1+5*cos(2*Ml)+el^3))/rol^2;fal2=2*sqrt(mu*rol*(1-el^2))*(el*nl*sin(Ml)+el^2*(nl*sin(2*Ml)+3*nl*cos(Ml)*sin(Ml))+el^3)/rol^2;c=diag([0.1;0.1;0.1]);% tao=diag([0.05;0.05;0.05]);beta=0.05;b=133*eye(3);% f0=[1;1;1];% D0=diag([0.5;0.5;0.5]);xp=[x(7);x(8);x(9)];xp1=[x(10);x(11);x(12)];xpd=[x(1);x(2);x(3)];xpd1=[x(4);x(5);x(6)];e=xp-xpd;e1=xp1-xpd1;s=e1+c*e;xpd2=[nl^2*x(1)+0*x(2)+2*nl*x(5)+2*nl^2*x(1); -0*x(1)+nl^2*x(2)-2*nl*x(4)-nl^2*x(2); -nl^2*x(3);];% xpd2=[0.1;-0.6;-0.5];XYZ=inv(A)*[rol+x(7);x(8);x(9)];X=XYZ(1,1);Y=XYZ(2,1);Z=XYZ(3,1);R=sqrt(X^2+Y^2+Z^2);%考虑地球J2摄动时随时间变化的位置半径gx1=-mu*x(7)/R^3*(1+1.5*J2*((Re/R)^2)*(1-5*(x(9)^2)/R^2));%地球J2摄动下X方向的地球引力项gy1=-mu*x(8)/R^3*(1+1.5*J2*((Re/R)^2)*(1-5*(x(9)^2)/R^2));%地球J2摄动下Y方向的地球引力项gz1=-mu*x(9)/R^3*(1+1.5*J2*((Re/R)^2)*(3-5*(x(9)^2)/R^2));%地球J2摄动下Z方向的地球引力项gxyz=A*[gx1;gy1;gz1];gx=gxyz(1,1);gy=gxyz(2,1);gz=gxyz(3,1);x0=[zeros(1,3),120];c1=[-0.5 0 0.5; -0.5 0 0.5];b=2;node=3;yd=[x(7) x(8) x(9)];dyd=[x(10) x(11) x(12)];ddyd=[2*nl*x(11)+3*nl^2*x(7) -2*nl*x(10) -nl^2*x(9)];e=[x(7)-x(1) x(8)-x(2) x(9)-x(3)];de=[x(10)-x(4) x(11)-x(5) x(12)-x(6)];kp=0.5;kd=0.3;K=[kp kd]';E=[e;de];Fai=[0 1;-kp -kd];A=Fai';Q=[500 0;0 500];P=lyap(A,Q);W=[x(10) x(11) x(12)]';xi=[e;de];h=zeros(3,1);for j=1:1:3 h(j)=exp(-norm(xi(:,j)-c1(:,j))^2/(2*b^2));end 四、运行结果
1 matlab版本
2014a
2 参考文献
[1] 门云阁.MATLAB物理计算与可视化[M].清华大学出版社，2013.