本文为控制数学的matlab求解这本书的笔记, 主要讲解控制问题的数学求解

大致上, 内容与数值分析等书一致, 主要的数学内容分为五大块, 分别为积分, 微分, 差分, 线代, 优化, 之后这本书依据这五个基本内容讨论了智能优化和鲁棒优化, 其中鲁棒优化主要是基于Riccati方程的求解上

第一波笔记

这是第一天的笔记, 打算分天把这本书看完.

% 将一个矩阵数值化
sym([[1,2],[2,3]]);
%%
%微积分与矩阵微积分运算
%二元微分函数的绘制
syms x y;
z = (x^2-x*2)*exp(-x^2-y^2-x*y);
dx = simplify(diff(z,x));%偏微分
dy = diff(z,y);%偏微分
latex(dx)%可以用latex函数将其转换成latex形式
[x1,y1] = meshgrid(-3:.2:3,-2:.2:2);
z1 = subs(z,{x,y},{x1,y1});
dx1 = subs(dx,{x,y},{x1,y1});
dy1 = subs(dy,{x,y},{x1,y1});
dx1 = double(dx1);
dy1 = double(dy1);
z1 = double(z1);%需要double()将sym()矩阵转换过来
subplot(1,3,1),surf(x1,y1,z1)
subplot(1,3,2),quiver(x1,y1,dx1,dy1)
subplot(1,3,3),contour(x1,y1,z1,30)
%%
%多重微分
syms x
f = exp(-5*x)*sin(2*x)/(x^2+3*x+5);
n = 7;
scope = -5:.02:-3;
fs = double(subs(f,x,scope));
plot(scope,fs)
hold on
for i = 1:n
    dx = diff(f,x,i);
    dxx = double(subs(dx,x,scope));
    plot(scope,dxx)
    hold on
end
%%
%隐函数的偏导, f(x,y), 求y对x的偏导, 有如下公式
sysm x y;
f = x^3+y^2*sin(x)-y-6;
dx = diff(f,x);
dy = diff(f,y);
dydx = -dx/dy;
%%
%参数方程多重导数
function result = paradiff(y,x,t,n)
if mod(n,1)~=0, error('Error: n should be positive integer');
else
    if n == 1
        result = diff(y,t)/diff(x,t);
    else
        result = diff(paradiff(y,x,t,n-1),t)/diff(x,t);
    end
end
end
%%
%LAPLACE变换, 用来将微分方程变换为多项式方程
syms t s t1;
f = t*exp(-3*t)*sin(2*t);
f1 = laplace(f)
%pretty(f1)
ilaplace(f1)
ilaplace(f1,s,t1)%符号用来转换不同的符号, 默认是复数s,时域t
%%
%对一个控制系统
%有tf, zpk, ss三种表达方式, 转换可以用tf2zp,zp2tf,ss2zp函数等等
%注意的是ioDelay可以描述传递函数乘以一个exp(-time)的函数
%egg
sys = tf(0.97,[1,2,34]);
sys1 = tf(0.97,[1,2,34],'ioDelay',0.73);
step(sys)
hold on
step(sys1)
%可以看到系统延时0.73s才进行响应, 原因是传递函数是在复数域上, 需要用exp来表述
%%
%对一个典型的二阶系统
%系统可以分为四种情况
%无阻尼, 欠阻尼, 临界阻尼, 过阻尼
%每种情况都有自己的阶跃响应表达式
wn = 1.5;
zeta = 0:.1:2;
t = 0:.05:20;
y = [];
for i = zeta
    [n,d] = ord2(wn,i);
    sys = tf(n,d);
    temp = step(sys,t)';
    if max(temp) >= 1
        flag = i;
    end
    y = [y;temp];
end
subplot(1,2,1),plot(t,y)
subplot(1,2,2),surf(t,zeta,y)
flag
%这里surf的用法很trick, 要记住这个用法, 在变量变化时对应不同映射的反映上很明显
%%
%eig与svd,
a = rand(3);
[V,D] = eig(a);%特征值分解, 使得a*V=V*D;
[U,S,V] = svd(a);%使得a=U*S*V';
%%
%提取矩阵特征多项式可以用poly, 只有方阵才有特征多项式, f(x) = det(sI-A)
a = rand(3);
pol = poly(a)
%求值用polyval与polyvalm(矩阵使用)
%vander矩阵, hilb矩阵, hankel矩阵, poly2sym函数,将p转换成函数形式.
A = vander([1 2 3 4 5 6 7]);
p = poly(A);
B = polyvalm(p,A)
norm(B)
%%
%秩, 迹
rank(a)
trace(a)
length(a)
size(a)

%%
%矩阵的逆
a = hilb(4);
inv1 = inv(a);
ans = inv1*a
err = eye(4)-ans
norm(err)%存在误差, 因此对于接近奇异的矩阵, 一般用pinv, 即广义逆矩阵, 也可以先把矩阵sym化, 然后求.
%广义逆对奇异矩阵也可行, 同时对非方阵也行
%定义为ANA=A, 其中N就是广义逆阵
%但是第一条定义很多矩阵都满足
%因此需要二三定义
%MAM=M
%AM与MA都是hermite对阵矩阵
%其定义为moore-penrose广义逆
%所谓hermite对称矩阵, 需要满足AM=M'A'
%%
%kronecker积
a = rand(3)
b = rand(4)
kron(a,b)
%%
%对一个二次型函数
%f(x)=x'*A*x-2*b'*x
%求一个导, 让导数为0即可.
%其极值可以转换成Ax=0形式, 直接左除即可
%%
%jacobian矩阵
%对有n个变量的m个函数
%可以用jacobian矩阵进行求导
syms r theta phi;
x = r*sin(theta)*cos(phi);
y = r*sin(theta)*sin(phi);
z = r*cos(theta);
J = jacobian([x y z],[r theta phi]);
jacobian(x,[r theta phi])
J
%可以想象成y行,x列的函数, y作为行需要作为列向量输入. 
%不过其实都按照行向量输入都是可以的
%不过没有hessian矩阵, 所以可以直接用两此jacobian来进行求解
%比如下
jacobian(jacobian(x,[r theta phi]),[r theta phi])