登录后复制
%%clc;clear;close all;warning off;N = 8;%%syms x1 x2;%initialize u,p,w 这个initial control law本来应该是给的,但是这个作者没给,随便设的我 Cost function不会有太大影响吧? u0 =-x1; %initialize state space,f,g,lx =[x1,x2];f =[-x1^3-x2;x1+x2];g =[0,140]';%相当于变更U0初始值l = x1'*x2;R = eye(1); %感觉应该是单位矩阵Ls = 1.5;%set the basis function phi and omega[-1,1]phi1 = x1^2; phi2 = x2*x1;phi3 = x2^2;phi4 = x1^4;phi5 = x1^3*x2;phi6 = x1^2*x2^2;phi7 = x1*x2^3;phi8 = x2^4;phi = [phi1,phi2,phi3,phi4,phi5,phi6,phi7,phi8];% N=8phi_diff = jacobian(phi,[x1 x2]); %partial differential of phi%%A_0 = int([phi_diff*f]'*phi',x2,-Ls,Ls)+int([phi_diff*g*u0]'*phi',x2,-Ls,Ls); b_0 =-int(l*phi,x2,-Ls,Ls)-int(u0'*R*u0*phi,x2,-Ls,Ls);c_0 = A_0^-1*b_0; %迭代部分syms c_1 c_2 c_3 c_4 c_5c{1} = c_0;c{2} = c_1;c{3} = c_2;c{4} = c_3;c{5} = c_4;M = 0;%外部迭代次数过多的话,会使得C的表达式越来越复杂,而导致仿真及其缓慢,这里根据我的电脑配置,设置了迭代3次for i=1:4 i Ctmp = c{i}; for j=1:N h = g*R^-1*g'; M = M + Ctmp(j)*int([phi_diff*h*[phi_diff(j,:)]']'*phi',x2,-Ls,Ls); end A = int([phi_diff*f]'*phi',x2,-Ls,Ls) - M/2; b =-int(l*phi,x2,-Ls,Ls)-(M/4)*c{i}; %这个c(i-1)是矩阵 c{i+1} = A^-1*b;end%obtain the cNC_N = A^-1*b;disp('u');syms UNUN = int(-1/2 * R * g' * phi_diff' * C_N',x2,-Ls,Ls);%可能是部分参数的问题,这里考虑到前面参数phi有4次的乘法,感觉需要最后开2次根号来抵消,所以做了0.25次方运算,%这样结果正好和论文一模一样了,所以这么做应该是合理的。UN = UN^0.25; index = 0;SCALE = [-1:0.1:1];y = [];for x1 = SCALE; x1 x1 = x1; index = index + 1; y(index) = abs(double(sum(subs(UN))));endfigure;plot(SCALE,y);xlabel('x_1');ylabel('Cost');save RN_8.mat SCALE y1.2.3.4.5.6.7.8.9.10.11.12.13.14.15.16.17.18.19.20.21.22.23.24.25.26.27.28.29.30.31.32.33.34.35.36.37.38.39.40.41.42.43.44.45.46.47.48.49.50.51.52.53.54.55.56.57.58.59.60.61.62.63.64.65.66.67.68.69.70.71.72.73.74.75.76.77.78.79.80.81.82.83.84.85.86.87.88.89.
免责声明:本文系网络转载或改编,未找到原创作者,版权归原作者所有。如涉及版权,请联系删