matlab of control engineers

一夕ξ

已于 2023-04-01 09:03:52 修改

阅读量1.8k

点赞数 1

分类专栏： matlab 文章标签： matlab 开发语言

于 2021-08-12 21:45:35 首次发布

本文链接：https://blog.csdn.net/qq_45387575/article/details/119650682

版权

matlab 专栏收录该内容

16 篇文章 1 订阅

订阅专栏

各大符号在matlab中的表示：

t=0:0.01*pi:2*pi
alpha=3
y=sin(alpha*t)
plot(t,y)
grid
title('plot of sin(\alphat)')
xlabel('t')
ylabel('sin(\alphat)')

u1=[0:0.5:2]'
u2=[1.5:-0.5:-2]'
u3=[-1.5:0.5:0]'
w=[u1;u2;u3]
t=0:0.5:8.0
plot(t,w)
grid

绘制两个矩形：

x1=[-10 10 10 -10 -10]
y1=[-5 -5 5 5 -5]
x2=[-5 5 5 -5 -5]
y2=[-2 -2 2 2 -2]
plot(x1,y1,x2,y2)

v=[-20 20 -20 20]
axis(v)

axis('square')

grid
xlabel('x')
ylabel('y')

xp=[0 5*sin(pi/6)]
yp=[0 5*cos(pi/6)]
plot(xp,yp)

v=[-6 6 -4 8]
axis(v)

hold
a=0:0.01:2*pi
x=5*sin(pi/6)+0.5*cos(a)
y=5*cos(pi/6)+0.5*sin(a)
plot(x,y)

xb=[-3 3 3 -3 -3]
yb=[-3 -3 0 0 -3]
plot(xb,yb)

plot(2+0.5*cos(a),-3.5+0.5*sin(a))

axis('square')

[r,p,k]=residue(num,den)：将传递函数进行部分分式化。

$\frac{2s^3+5s^2+3s+6}{s^3+6s^2+11s+6}=\frac{-6}{s+3}+\frac{-4}{s+2}+\frac{3}{s+1}+2$

num=[2 5 3 6]
den=[1 6 11 6]
[r,p,k]=residue(num,den)

>> main

num =

2 5 3 6

den =

1 6 11 6

r =

-6.0000
-4.0000
3.0000

p =

-3.0000
-2.0000
-1.0000

k =

2

也可以反过来使用：

[num,den]=residue(r,p,k)

r=[-6 -4 3]
p=[-3 -2 -1]
k=2
[num,den]=residue(r,p,k)

>>

num =

2 5 3 6

den =

1 6 11 6

$\frac{B(s)}{A(s)}=\frac{num}{den}$

>> printsys(num,den,'s')

num/den =

2 s^3 + 5 s^2 + 3 s + 6
-----------------------
s^3 + 6 s^2 + 11 s + 6

num=[1 2 3]
den=[1 3 3 1]
[r,p,k]=residue(num,den)

>>

r =

1.0000
0.0000
2.0000

p =

-1.0000
-1.0000
-1.0000

k =

[]

$\frac{s^2+2s+3}{s^3+3s^2+3s+1}=\frac{1}{s+1}+\frac{0}{(s+1)^2}+\frac{2}{(s+1)^3}$

[z,p,k]=tf2zp(num,den):求零极点表达式

z是零点，p是极点，k是增益

num=[4 16 12]
den=[1 12 44 48]
[z,p,k]=tf2zp(num,den)

>>

z =

-3
-1

p =

-6.0000
-4.0000
-2.0000

k =

4

[num,den]=zp2tf(z,p,k);printsys(num,den,'s'):将零极点表达换为一般表达;z,p一定是列向量

z=[-1;-3]
p=[0;-2;-4;-6]
k=4
[num,den]=zp2tf(z,p,k)
printsys(num,den,'s')

>>

num =

0 0 4 16 12

den =

1 12 44 48 0

num/den =

4 s^2 + 16 s + 12
----------------------------
s^4 + 12 s^3 + 44 s^2 + 48 s

num=[5 330 55 30]
den=[1 9 33 65]
[z,p,k]=tf2zp(num,den)

>>

z =

-65.8343 + 0.0000i
-0.0829 + 0.2903i
-0.0829 - 0.2903i

p =

-5.0000 + 0.0000i
-2.0000 + 3.0000i
-2.0000 - 3.0000i

k =

5
$\frac{2s^3+330s^2+55s+30}{s^3+9s^2+33s+65}=\frac{5(s+65.8343)(s+0.0829-0.2903i)(s+0.0829+0.2903i)}{(s+5)(s+2-3i)(s+2+3i)}$

tf2ss

ss2tf

ss2zp

zp2ss

tf2zp

zp2tf

c2d

num=[25.04 5.008]
den=[1 5.03247 25.1026 5.008]
[A,B,C,D]=tf2ss(num,den)
>>

A =

-5.0325 -25.1026 -5.0080
1.0000 0 0
0 1.0000 0

B =

1
0
0

C =

0 25.0400 5.0080

D =

0

[num,den]=ss2tf(A,B,C,D,iu):对于单输出多输入系统，由状态方程转换成传递函数来说有多个

A=[0 1;-2 -3]
B=[1 0;0 1]
C=[1 0]
D=[0 0]
[num1,den1]=ss2tf(A,B,C,D,1)
[num2,den2]=ss2tf(A,B,C,D,2)

>>

>> main

A =

0 1
-2 -3

B =

1 0
0 1

C =

1 0

D =

0 0

num1 =

0 1 3

den1 =

1 3 2

num2 =

0 0 1

$x'=\begin{pmatrix} 0&1 \\ -2&-3 \end{pmatrix}x+\begin{pmatrix} 1&0 \\ 0&0 \end{pmatrix}u$

$y=\bigl(\begin{smallmatrix} 1 & 0 \end{smallmatrix}\bigr)x$

$\frac{U}{R1}=\frac{s+3}{s^2+3s+2}$

$\frac{U}{R2}=\frac{1}{s^2+3s+2}$

[num,den]=ss2tf(A,B,C,D,iu):对于输出多输入系统

A=[0 1;-25 -4]
B=[1 1;0 1]
C=[1 0;0 1]
D=[0 0;0 0]
[num1,den1]=ss2tf(A,B,C,D,1)
[num2,den2]=ss2tf(A,B,C,D,2)

num1 =

0 1 4
0 0 -25

den1 =

1.0000 4.0000 25.0000

num2 =

0 1.0000 5.0000
0 1.0000 -25.0000

den2 =

1.0000 4.0000 25.0000

$x'=\bigl(\begin{smallmatrix} 0 & 1\\ -25& -4 \end{smallmatrix}\bigr)x+\bigl(\begin{smallmatrix} 1 & 1\\ 0& 1 \end{smallmatrix}\bigr)u$

$y=\bigl(\begin{smallmatrix} 1 &0 \\ 0 & 1 \end{smallmatrix}\bigr)x+\bigl(\begin{smallmatrix} 0 &0 \\ 0& 0 \end{smallmatrix}\bigr)u$

$\frac{U1}{R1}=\frac{s+4}{s^2+4s+25}$

$\frac{U2}{R1}=\frac{-25}{s^2+4s+25}$

$\frac{U2}{R1}=\frac{s+5}{s^2+4s+25}$

$\frac{U2}{R2}=\frac{s-25}{s^2+4s+25}$

minreal(sys):将表达式sys中的零极点相消后再显示

num=[1 3 2]
den=[1 8 19 12]
G=minreal(tf(num,den))

>>

G =

s + 2
--------------
s^2 + 7 s + 12

响应：

wn=5
zeta=0.4
num=[wn^2]
den=[1 2*zeta*wn wn^2]
G=tf(num,den)

>>

G =

25
--------------
s^2 + 4 s + 25

单位阶跃响应：

step(sys,t);step(num,den,t);y=step(sys,t);y=step(num,den,t)[y,t]=step(sys,t);[y,t]=step(num,den,,t)

t=0:0.1:10
[y,t]=step(G,t)
plot(t,y)
grid

m=10,b=20,k=100

num=[2 10]
den=[1 2 10]
t=0:0.1:10
[y,t]=step(num,den,t)
plot(y)
grid

>>

>> printsys(num,den,'s')

num/den =

2 s + 10
--------------
s^2 + 2 s + 10

axis([0 60 0 1.8])
title('step response')
xlabel('t(sec)')
ylabel('y')

xlabel('string','Fontsize',20):

axis([0 60 0 1.8])
title('step response','Fontsize',20)
xlabel('t(sec)')
ylabel('y')

多输入多输出系统的单位阶跃响应

A=[-1 -1;6.5 0]
B=[1 1;1 0]
C=[1 0;0 1]
D=[0 0;0 0]
step(A,B,C,D)

将同一个输入下的单位阶跃响应图绘制再一张图上：

step(A,B,C,D,1)
gtext('in1-out1')
gtext('in1-out2')

t=0:0.01:5
num1=[1]
den1=[1 2 0]
sys1=tf(num1,den1)
sys=feedback(10*sys1,1)
[y,t]=step(sys,t)
plot(t,y)
grid

绘制单位阶跃响应图：

z=[]
p=[-1 -4]
k=20
sys1=zpk(z,p,k)
sys2=feedback(sys1,0.4490,-1)

num=[1]
den=[1 0]
sys3=tf(num,den)
sys4=series(sys2,sys3)

sys4=feedback(sys4,1,-1)

>>

sys4 =

20
--------------------------------
(s+2.902) (s^2 + 2.098s + 6.892)

t=0:0.01:8
[y,t]=step(sys4,t)
plot(t,y)
grid

feedback(sys1,sys2,-1/+1):表示正负反馈

feedback(sys1,1):表示单位负反馈

s=tf('s')
sys1=2/s
sys2=2/(s+4)
sys3=series(sys1,sys2)
sys4=feedback(sys3,1,-1)
sys5=0.75/s
sys6=series(sys5,sys4)

sys7=feedback(sys6,1,-1)

t=0:0.01:10
[y,t]=step(sys7,t)
plot(t,y)
title('step response','Fontsize',20)
xlabel('t(sec)')
ylabel('c(t)')

t=0:0.01:20
k=[5 10 20]
for n=1:3
numg=[k(n)];deng=[1 11 10 0]
sysg=tf(numg,deng)
sys=feedback(sysg,[1])
y=step(sys,t)
plot(t,y)
hold on
end

t=0:0.01:20
k=[5 10 20]
z=[]
p=[0 -1 -10]
for i=1:3
k2=[k(i)] 矩阵元素的访问
sys1=zpk(z,p,k2)
sys2=feedback(sys1,1,-1)
y(:,i)=step(sys2,t) 注意写法，数组的索引必须是正整数或者逻辑值
plot(t,y(:,i))
hold on
end
gtext('k=5')
gtext('k=10')
gtext('k=20')

subplot(2,2,1)
plot(t,y(:,1)),grid,title('k=5')
subplot(2,2,2)
plot(t,y(:,2)),grid,title('k=10')
subplot(2,2,3)
plot(t,y(:,3)),grid,title('k=20')

t=0:0.01:20
k=[0.1 0.2 0.3 0.4 0.5]
for i=1:5
s=tf('s')
sys1=12.5/(s+1)
k1=[k(i)]
sys2=feedback(sys1,k1,-1)
sys3=1/s
sys4=series(sys2,sys3)
sys5=feedback(sys4,1,-1)
y(:,i)=step(sys5,t)
plot(t,y(:,i))
hold on
end
gtext('k=0.1')
gtext('k=0.2')
gtext('k=0.3')
gtext('k=0.4')
gtext('k=0.5')

t=0:0.01:20
k=[0.1 0.2 0.3 0.4 0.5]
for i=1:5
s=tf('s')
sys1=12.5/(s+1)
k1=[k(i)]
sys2=feedback(sys1,k1,-1)
sys3=1/s
sys4=series(sys2,sys3)
sys5=feedback(sys4,1,-1)
y(:,i)=step(sys5,t)
end
plot(t,y)
gtext('k=0.1')
gtext('k=0.2')
gtext('k=0.3')
gtext('k=0.4')
gtext('k=0.5')

？？？？？？？？不会系统的近似

阶跃响应

状态空间中的脉冲响应，[y,x,t]:会有x的向量产生

A=[0 1;-1 -1]
B=[0;1]
C=[1 0]
D=[0]
t=0:0.01:10
[y,x,t]=impulse(A,B,C,D)
subplot(2,2,1)
plot(t,x(:,1))
subplot(2,2,2)
plot(t,x(:,2))
subplot(2,2,3)
plot(t,y)

num=[1]
den=[1 0.2 1]
t=0:0.01:10
y=impulse(num,den)
plot(y)

t=0:0.01:10

[t,y]=impulse(num,den)
plot(t,y)

>>

错误使用 plot
向量长度必须相同。

出错 main (第 5 行)
plot(t,y)

很奇怪，明明应该向量的长度就是一样的。

斜坡信号：ramp:需要自己设，用lsim指令：

对系统： $\frac{1}{s^2+s+1}$ 求其单位斜坡响应相当于对系统 $\frac{1}{s^2+s+1}\frac{1}{s}$ 求一个单位阶跃响应

num=[1]
den=[1 1 1 0]
t=0:0.01:7
c=step(num,den,t)

plot(t,c,'--',t,t,'g')

num1=[10 4]
den1=[1 4 4]
t=0:0.01:10
y1=step(num1,den1,t)
subplot(1,2,1)
plot(y1)
grid
title('step-response')

num2=[10 4]
den2=[1 4 4 0]
y2=step(num2,den2,t)
subplot(1,2,2)
plot(y2)
grid
title('ramp-response')

plot(t,y)

lsim(sys,u,t,x0):x0:为初始条件

num=[1]
den=[1 1 1]
sys=tf(num,den)
t=0:0.01:8
u=t 此时u为单位斜坡信号
[y,t]=lsim(sys,u,t)
plot(t,y)

A=[-1 0.5;-1 0]
B=[0;1]
C=[1 0]
D=[0]
sys=ss(A,B,C,D)
t=0:0.01:12
u1=t
u2=exp(-t)
x0=[0 0]
[y1,t]=lsim(sys,u1,t,x0)
[y2,t]=lsim(sys,u2,t,x0)
y3=step(sys,t)
plot(t,y1,'--',t,y2,'g--',t,y3,'r--')
axis([0 12 0 1.4])

num=[1 10]
den=[1 6 9 10]
sys=tf(num,den)
t=0:0.01:10
u1=t
u2=exp(-0.5.*t)
[y1,t]=lsim(sys,u1,t)
[y2,t]=lsim(sys,u2,t)
plot(t,y1,'--',t,y2,'g--')

legend('ramp-response','exp(-0.5t)-reaponse')

num=[5]
den=[1 1 5]
sys=tf(num,den)

t=0:0.01:8
r=2+t
[y,t]=lsim(sys,r,t)
plot(t,y)

num=[2]
den=[1 1]
sys=tf(num,den)
sys1=feedback(sys,1,-1)
t=0:0.01:10
u=0.5*t.^2
[y,t]=lsim(sys1,u,t)
plot(t,y)
xlabel('t(sec)')
ylabel('output')
title('input=0.5t^2--response','Fontsize',12)

A=[0 1;-5 -3]
B=[0;1]
C=[1 0]
D=0
sys=tf(ss(A,B,C,D))
t=0:0.01:10
u=0.5*t.^2
x0=[0;0]
[y,x]=lsim(sys,u,t,x0) 吧
plot(t,y)

练习：

？？？？？？？matlab中怎么编程：

A=[0 1;-10 -5]
B=[2;1]
t=0:0.01:10
[x,z,t]=step(A,B,A,B,1,t)
subplot(1,2,1)
plot(t,x(:,1))
subplot(1,2,2)
plot(t,x(:,2))

也可以使用[y,x]=initial(A,B,C,D,[initial condition],t)表示零状态响应

A=[0 1;-10 -5]
t=0:0.01:8
B=[0;0]
C=[0 0]
D=0
x0=[2;1]
[y,x]=initial(A,B,C,D,x0,t)
subplot(1,2,1)
plot(t,x(:,1))
subplot(1,2,2)
plot(t,x(:,2))

A=[0 1 0;0 0 1;-10 -17 -8]
B=[0;0;0]
C=[1 0 0]
D=0
t=0:0.01:10
x0=[2;1;0.5]
[y,x]=initial(A,B,C,D,x0,t)
subplot(2,2,1)
plot(t,x(:,1))
subplot(2,2,2)
plot(t,x(:,2))
subplot(2,2,3)
plot(t,x(:,3))

subplot(2,2,4)
plot(t,y)

zeta=[0 0.2 0.4 0.6 0.8]
t=0:0.01:10
for i=1:5
num=[1]
den=[1 2*zeta(i) 1]
sys=tf(num,den)
y=step(sys,t)
plot(t,y)
grid
hold on

end
legend('zeta=0')
hold on
legend('zeta=0.2')
legend('zeta=0.4')
legend('zeta=0.6')
legend('zeta=0.8')

????????不用gtext，怎么把每个曲线的标签加进去

zeta=[0 0.2 0.4 0.6 0.8]
t=0:0.01:10
for i=1:5
num=[1]
den=[1 2*zeta(i) 1]
sys=tf(num,den)
y(:,i)=step(sys,t)
end
plot(t,y)

mesh(zeta,t,y)

xlabel('\zeta')
ylabel('t(sec)')
zlabel('output')

mesh(t,zeta,y')
xlabel('t')
ylabel('\zeta')
zlabel('output') 输出都一样，只是此时需要把x,y坐标轴对换

root-locus analysis

r=locus(num,den)

plot(r,'--')

num=[1 3]

a=[1 1 0]
b=[1 4 16]
c=conv(a,b)
den=[1 5 20 16 0]
sys=tf(num,den)
r=rlocus(num,den)
v=[-6 6 -6 6]
pzmap(sys)
hold on
plot(r,'r-')
grid
axis(v)

num=[1]
den=conv([1 0],conv([1 0.5],[1 0.6 10]))
sys=tf(num,den)
pzmap(num,den)
hold on
r=rlocus(num,den)
plot(r,'r-')
v=[-6 6 -6 6]
axis(v)

num=[1]
den=conv([1 0],conv([1 0.5],[1 0.6 10]))
sys=tf(num,den)
pzmap(num,den)
hold on
k1=0:0.2:20
r1=rlocus(num,den,k1)
plot(r1,'r-')
k2=[20:0.1:30]
k3=[30:0.5:1000]
r2=rlocus(num,den,k2)
r3=rlocus(num,den,k3)
plot(r2,'g--')
plot(r3,'b-')
v=[-6 6 -6 6]
axis(v)

A=[0 1 0;0 0 1;-160 -56 -14]
B=[0;1;-14]
C=[1 0 0]
D=[0]
k=0:0.1:400
r=rlocus(A,B,C,D,k)
plot(r,'g-')
hold on
pzmap(A,B,C,D)
v=[-16 1 -8 8]
axis(v)

asymptotes：渐近线

渐近线是s-∞大时的根轨迹

$\lim_{s\rightarrow\infty }\frac{k}{s^3+3s^2+2s}=\lim_{s\rightarrow\infty }\frac{k}{s^3+3s^2+3s+1}=\frac{k}{(s+1)^3}$

num=[1]
den=conv([1 1],conv([1 1],[1 1]))
r=rlocus(num,den)
plot(r,'b--')
hold on

num1=[1]
den1=conv([1 0],conv([1 1],[1 2]))
pzmap(num1,den1)
hold on
r2=rlocus(num1,den1)
plot(r2,'r-')
v=[-8 8 -8 8]
axis(v)

num=[1]
den=conv([1 1],conv([1 1],[1 1]))
r=rlocus(num,den)

num1=[1]
den1=conv([1 0],conv([1 1],[1 2]))
r2=rlocus(num1,den1)
plot([r r2])
v=[-8 8 -8 8]
axis(v)

错误使用 horzcat
要串联的数组的维度不一致。

>>为了保证r和r2矩阵的维度一致，需要增加k的变化

num=[1]
den=conv([1 1],conv([1 1],[1 1]))
k=[0:0.01:100]
r=rlocus(num,den,k)

num1=[1]
den1=conv([1 0],conv([1 1],[1 2]))
r2=rlocus(num1,den1,k)
plot([r r2])
v=[-8 8 -8 8]
axis(v)

num=[1]
den=conv([1 1],conv([1 1],[1 1]))
k1=[0:0.01:100]
k2=[100:0.2:400]
k3=[400:0.5:600]
k=[k1 k2 k3]
r=rlocus(num,den,k)

num1=[1]
den1=conv([1 0],conv([1 1],[1 2]))
pzmap(num1,den1)
hold on
r2=rlocus(num1,den1,k)
plot([r r2])
v=[-8 8 -8 8]
axis(v)

num=conv([1 2],[1 2])
den=conv([1 0 4],conv([1 5],[1 5]))
r=rlocus(num,den)
pzmap(num,den)
hold on
plot(r,'g-')
axis([-8 8 -8 8])

sgrid([a b...],[c d ...]):绘图，显示阻尼比和无阻尼振荡频率

sgrid([0.5,0.707],[0.5 1 2])
v=[-3 3 -3 3]
axis(v)
axis('square')

num=[0 0 0 1];
den=[1 4 5 0]
pzmap(num,den)
hold on
r=rlocus(num,den)
plot(r,'g-')
sgrid([0.5 0.7],[0.5 1 2])
v=[-3 3 -3 3]
axis(v)

sgrid(0.5,[])：只显示一部分

num=[1 2]
den=[1 9 8 0]
pzmap(num,den)
hold on
r=rlocus(num,den)
plot(r,'r-')
v=[-10 10 -10 10]
axis(v)
sgrid
title('r-locus')

[k,r]=rlocfind(num,den)或者[k,r]=rclofind(A,B,C,D)：寻找根轨迹中某一个特定点处的k值、以及其他的极点。

[k,r]=rlocfind(num,den)

>> [k,r]=rlocfind(num,den)
Select a point in the graphics window

selected_point =

-4.6023 - 0.0496i

k =

21.6480

r =

-4.6018 + 0.0000i
-2.1991 + 2.1383i
-2.1991 - 2.1383i

num=[1]
den=[1 4 5 0]
pzmap(num,den)
hold on
r=rlocus(num,den)
plot(r,'g-')
v=[-5 5 -5 5]
axis(v)
sgrid
[k,r]=rlocfind(num,den)

>>

Select a point in the graphics window

selected_point =

-0.7292 + 0.8685i

k =

3.2739

r =

-2.5383 + 0.0000i
-0.7308 + 0.8693i
-0.7308 - 0.8693i

绘制根轨迹是根据开环传递函数来的，而响应是对于某一个系统即闭环传递函数来的。

num1=[1 2]
den1=[1 4 0]
subplot(2,2,1)
pzmap(num1,den1)
hold on
r1=rlocus(num1,den1)
plot(r1,'g-')

num=[-1 2]
den=[1 4 0]
subplot(2,2,2)
pzmap(num,den)
hold on
r2=rlocus(num,den)
plot(r2,'r-')
v=[-5 5 -5 5]
axis(v)

sys=tf(num1,den1)
sys1=feedback(sys,1,-1)
sys2=tf(num,den)
sys3=feedback(sys2,1,-1)
t=0:0.01:10
y1=step(sys1,t)
y2=step(sys3,t)
subplot(2,2,3)
plot(t,y1)

subplot(2,2,4)
plot(t,y2)

绘制参数在特定变化范围内的根轨迹：

num=[1 2 4]
den=conv([1 0],conv([1 4],conv([1 6],[1 1.4 1])))
k1=[0:0.001:12]
k2=[73:0.5:154]
k=[k1 k2]
pzmap(num,den)
hold on
r=rlocus(num,den,k)
plot(r,'r-')
v=[-5 5 -5 5]
axis(v)
grid

pade近似：可以用一个多项式的传递函数对特定周期的延迟时间进行近似

pade(a,b):a表示延迟时间,b表示近似阶数

T=0.1
[num,den]=pade(0.1,3)
sys=tf(num,den)

>>

sys =

-s^3 + 120 s^2 - 6000 s + 1.2e05
--------------------------------
s^3 + 120 s^2 + 6000 s + 1.2e05

Continuous-time transfer function.

$e^{-0.1s}=\frac{-s^3+120s^2-6000s+120000}{s^3+120s^2+6000s+120000}$

specification：规定、规范

requency-response anaylysis

mag就是普通的单位，要绘制伯德图的话，需要计算L(w)=20*log10(mag)

semilogx(w,20*log10(mag)):才能获得正确的伯德图

用w=logspace(a,b):产生w的段位10^a~~10^b之间的50个点，或者w=logspace(a,b,n):产生n个点

num=[25]
den=[1 4 25]
bode(num,den)

num=[25]
den=[1 4 25]
w=logspace(-1,3)
bode(num,den,w)

num=[9 1.8 9]
den=[1 1.2 9 0]

[mag,phase,w]=bode(num,den,w)
w=logspace(-2,3,100)
[mag,phase,w]=bode(num,den,w) 并不会有图像生成

magdB=20*log10(mag)
dBmax=50*ones(1,100)
dBmin=-50*ones(1,100)

semilogx(w,magdB,'o',w,magdB,'-',w,dBmax,'-',w,dBmin,'-')

grid
title('Bode diagram')
xlabel('frequency(rad/sec)')
ylabel('gain dB')

pmax=150*ones(1,100)
pmin=-150*ones(1,100)
semilogx(w,phase,'o',w,phase,'-',w,pmax,'-',w,pmin,'-')

bode(A,B,C,D):会产生多个伯德图

A=[0 1;-25 -4]
B=[0;25]
C=[1 0]
D=0
bode(A,B,C,D)/bode(A,B,C,D,1)

A=[0 1;-25 -4]
B=[1 1;0 1]
C=[1 0;0 1]
D=0
bode(A,B,C,D)

bode(A,B,C,D,iu):绘制第iu个输入对全部输出的bode图

bode(A,B,C,D,1)

k=[1 10 20]
for i=1:3
num=[k(i)]
den=conv([1 0],conv([1 1],[1 5]))
sys=tf(num,den)
sys1=feedback(sys,1,-1)
[mag,phase,w]=bode([k(i)],[1 6 5 k(i)],w)
subplot(1,2,1)
magdB=20*log10(mag)
semilogx(w,magdB)
hold on
subplot(1,2,2)
semilogx(w,phase)
hold on
end

行不通，mag始终是一个三维数组。？？？？？？？？？？

w=logspace(-1,2,200)
for i=1:3
if i==1;k=1;[mag,phase,w]=bode([k],[1 6 5 k],w);
mag1dB=20*log10(mag);phase1=phase;end
if i==2;k=2;[mag,phase,w]=bode([k],[1 6 5 k],w);
mag2dB=20*log10(mag);phase2=phase;end
if i==3;k=3;[mag,phase,w]=bode([k],[1 6 5 k],w);
mag3dB=20*log10(mag);phase3=phase;end
end
semilogx(w,mag1dB,'-',w,mag2dB,'-',w,mag3dB,'-')
grid
gtext('k=1')
gtext('k=2')
gtext('k=3')

semilogx(w,phase1,'-',w,phase2,'-',w,phase3,'-')
grid
gtext('k=1')
gtext('k=2')
gtext('k=3')

用subplot绘制：

奈奎斯特图：

此时只能生成re,im,w的矩阵，并不会有曲线生成

num=[1]
den=[1 0.8 1]
nyquist(num,den)

v=[-2 2 -2 2]
axis(v)
title('nyquist map')

合理改变坐标轴的范围，可以获得有用的信息

num=[1]
den=[1 1 0]
w=0.1:0.1:100
[re,im,w]=nyquist(num,den,w)
plot(re,im)
v=[-2 2 -2 2]
grid
axis(v)
title('nyquist map')

A=[0 1;-25 -4]
B=[0;25]
C=[1 0]
D=0
nyquist(A,B,C,D)

A=[-1 -1;6.5 0]
B=[1 1;1 0]
C=[1 0;0 1]
D=0
nyquist(A,B,C,D)

num=[20 20 20*0.5]
den=conv([1 0],conv([1 1],[1 10]))
w1=0.1:0.1:10
w2=10:0.1:100
w3=100:0.1:400
w=[w1 w2 w3]
[re,im,w]=nyquist(num,den,w)
plot(re,im)

找出特定点处对应的实部虚部、幅值、相角

num=[20 20 20*0.5]
den=conv([1 0],conv([1 1],[1 10]))
w1=0.1:0.1:10
w2=10:0.1:100
w3=100:0.1:400
w=[w1 w2 w3]
[re,im,w]=nyquist(num,den,w)
plot(re,im)
hold on
w1=[0.2 0.3 0.5 1 2 6 10 20]
[re,im,w]=nyquist(num,den,w1)
plot(re,im,'o')
v=[-2 2 -2 2]
axis(v)
grid

[re,im]
>>

>> [re,im]

ans =

0.9419 -4.8265
0.9899 -3.0878
1.1172 -1.6559
1.4356 -0.6436
1.7115 -0.4423
1.4487 -0.8737
0.9946 -0.9955
0.3995 -0.7990

[mag,phase,w]=bode(num,den,w1)
[mag,phase,w]

>>

ans =

4.9176 -78.9571 0.2000
3.2426 -72.2244 0.3000
1.9975 -55.9925 0.5000
1.5733 -24.1455 1.0000
1.7678 -14.4898 2.0000
1.6918 -31.0946 6.0000
1.4072 -45.0285 10.0000
0.8933 -63.4385 20.0000

k=[1 7.5 20]
for i=1:3
if i==1;k=1;num=[10*k];den=conv([1 0],conv([1 5],[0.1 1]));[re1,im1]=nyquist(num,den)
end
if i==2;k=7.5;num=[10*k];den=conv([1 0],conv([1 5],[0.1 1]));[re2,im2]=nyquist(num,den)
end
if i==3;k=20;num=[10*k];den=conv([1 0],conv([1 5],[0.1 1]));[re3,im3]=nyquist(num,den)
end
end
plot(re1,im1,'r-',re2,im2,'g-',re3,im3,'b-')
axis([-20 1 -20 20])
gtext('k=1')
gtext('k=7.5')
gtext('k=20')

？？？？？？？代码始终运行不完

num=[20 20]
den=conv([0 1 0],conv([0 1 5],[1 2 10]))
sys=tf(num,den)
bode(sys)
hold on
[Gm,Pm,wcp,wcg]=margin(sys)

>>

Gm =

3.1369

Pm =

103.6573

wcp =

4.0132

wcg =

0.4426

matlab approach to state-space design of control system

acker()：进行极点配置是基于艾克尔曼公式来的，其只能用于单输入的系统

place(A,B,p):缺点p的秩不能比B的秩大或者等于。

A=[0 1 0;0 0 1;-1 -5 -6]
B=[0;0;1]
p=[-2+4j -2-4j -10]
l=rank(ctrb(A,B))
k=place(A,B,p)
k2=acker(A,B,p)

>>

l =

3

k =

199.0000 55.0000 8.0000

k2 =

199 55 8

引入的状态反馈是u=-kx

A=[0 1 0;0 0 1;-1 -5 -6]
B=[0;0;1]
C=[1 0 0]
D=0
t=0:0.1:10
x0=[1;0;0]

[y,x]=initial(A,B,C,D,x0,t)
subplot(2,2,1)
plot(t,x(:,1))
subplot(2,2,2)
plot(t,x(:,2))
subplot(2,2,3)
plot(t,x(:,3))

A=[0 1 0;0 0 1;-1 -5 -6]
B=[0;0;1]
C=[1 0 0]
D=0
t=0:0.1:20
x0=[1;0;0]
k=[199 55 8]
A2=A-B*k
B2=[0;0;0]
[y,x]=initial(A,B,C,D,x0,t)
[y2,x2]=initial(A2,B2,C,D,x0,t)
subplot(2,2,1)
plot(t,x(:,1),'b-',t,x2(:,1),'r-')
subplot(2,2,2)
plot(t,x(:,2),'b-',t,x2(:,2),'r-')
subplot(2,2,3)
plot(t,x(:,3),'b-',t,x2(:,3),'r-')

矩阵C,D并不会影响状态反馈增益矩阵K

A=[0 1 0;0 0 1;-6 -11 -6]
B=[0;0;10]
p=[-2+2*sqrt(3)*j -2-2*sqrt(3)*j -10]
l=rank(ctrb(A,B))
k=place(A,B,p)
k2=acker(A,B,p)

>>

l =

3

k =

15.4000 4.5000 0.8000

k2 =

15.4000 4.5000 0.8000

以上极点配置时，都是假设完全可控的系统中所有的状态变量都是可以用来反馈的。但在实际上，并不是所有的状态变量都能用来反馈，因此我们需要对不能反馈的状态变量进行观测。设计状态观测器---

但对于输出我们总是能够观测到，所以与输出有关的状态变量不需要额外设计状态观测器，因此总的需要观测状态变量的个数位n-m，n为状态向量的维数，m为输出向量的维数。

全维观测器：

状态观测的两个关键点在于：‘复制’、‘反馈’

L=acker(A',C',P)' L=place(A',C',P)'

P是状态观测器的期望极点

降维观测器：

L=acker(Abb',Aab',P)' L=place(Abb',Aab',P)'

状态反馈增益矩阵：

A=[0 1 0;0 0 1;-6 -11 -6]
B=[0;0;1]
p=[-2+2*sqrt(3)*j -2-2*sqrt(3)*j -6]
l=rank(ctrb(A,B))
k=place(A,B,p)

>>

l =

3

k =

90.0000 29.0000 4.0000

设置观测器，由于y=cx,所以x1是可以观测到的，只需要观测x2,x3-设置降维观测器

Abb=[0 1;-11 -6]
A1=Abb'
Aab=[1 0]
A2=Aab'
p=[-10 -10]
k=acker(A1,A2,p)
L=k'

>>

k =

14 5

L =

14
5

A=[0 0 1 0;0 0 0 1;-36 36 0.6 -0.6;18 -18 0.3 -0.3]
B=[0;0;1;0]
C=[1 0 0 0;0 1 0 0]
D=0

l=rank(ctrb(A,B))
P=[-2+2*sqrt(3)*j -2-2*sqrt(3)*j -10 -10]
K=acker(A,B,P)

>>

K =

130.4444 -41.5556 24.3000 14.2185

Aab=[1 0;0 1]
Abb=[0.6 -0.6;0.3 -0.3]
P2=[-15 -16]
L=place(Abb',Aab',P2)'

>>

L =

15.6000 -0.6000
0.3000 15.7000

F=[0 0;0 0;1 0;0 1]
AA=[A-B*K B*K*F;zeros(2,4) Abb-L*Aab]
sys=ss(AA,eye(6),eye(6),eye(6))
t=0:0.1:4
x0=[0.1 0 0 0 0.1 0.05]
x=initial(sys,x0,t) [y,x]=initial(A,B,C,D,x0,t):要分清变量和输出曲线的对应？？？？？？？？？？？？？？？
subplot(2,3,1)
plot(t,x(:,1)),grid,xlabel('t'),ylabel('x1')
subplot(2,3,2)
plot(t,x(:,2)),grid,xlabel('t'),ylabel('x2')
subplot(2,3,3)
plot(t,x(:,3)),grid,xlabel('t'),ylabel('x3')
subplot(2,3,4)
plot(t,x(:,4)),grid,xlabel('t'),ylabel('x4')
subplot(2,3,5)
plot(t,x(:,5)),grid,xlabel('t'),ylabel('e1')
subplot(2,3,6)
plot(t,x(:,6)),grid,xlabel('t'),ylabel('e2')

未进行状态反馈时的零状态响应曲线：

A=[0 0 1 0;0 0 0 1;-36 36 0.6 -0.6;18 -18 0.3 -0.3]
B=[0;0;1;0]
C=[1 0 0 0;0 1 0 0]
D=[0;0]
x0=[0.1 0 0 0]
t=0:0.1:10
[y,x,t]=initial(A,B,C,D,x0,t)
subplot(2,2,1)
plot(t,x(:,1)),grid,xlabel('t'),ylabel('x1')
subplot(2,2,2)
plot(t,x(:,2)),grid,xlabel('t'),ylabel('x2')
subplot(2,2,3)
plot(t,x(:,3)),grid,xlabel('t'),ylabel('x3')
subplot(2,2,4)
plot(t,x(:,4)),grid,xlabel('t'),ylabel('x4')

A=[0 1;20.6 0]
B=[0;1]
C=[1 0]
P=[-1.8+2.4j -1.8-2.4j]
K=place(A,B,P)

>>

K =

29.6000 3.6000

P2=[-8 -8]
L=acker(A',C',P2)'

>>

L =

16.0000
84.6000

L就是Ke

num=[1]
den=[1 0 -20.6]
sys1=tf(num,den)

num1=[778.2 3690.7]
den1=[1 19.6 151.2]
sys2=tf(num1,den1)
sys3=feedback(sys1,sys2,-1)

>>
$\frac{Y(s)}{U(s)}$

sys3 =

s^2 + 19.6 s + 151.2
------------------------------------------
s^4 + 19.6 s^3 + 130.6 s^2 + 374.4 s + 576

A=[0 1;20.6 0]
K =[29.6000 3.6000]
L=[16;84.6]
B=[0;1]
C=[1 0]
AA=[A-B*K B*K;zeros(2,2) A-L*C]
x0=[1;0;0.5;0]
t=0:0.01:10
x=initial(AA,eye(4),eye(4),eye(4),x0,t)
subplot(2,2,1)
plot(t,x(:,1)),grid,xlabel('t'),ylabel('x1')
subplot(2,2,2)
plot(t,x(:,2)),grid,xlabel('t'),ylabel('x2')
subplot(2,2,3)
plot(t,x(:,3)),grid,xlabel('t'),ylabel('e1')
subplot(2,2,4)
plot(t,x(:,4)),grid,xlabel('t'),ylabel('e2')

基于降维观测器的状态反馈设计：

解题步骤：

A=[0 1 0;0 0 1;0 -24 -10]
B=[0;10;-80]
C=[1 0 0]
D=0
[num,den]=ss2tf(A,B,C,D)

>>

num =

0 0 10.0000 20.0000

den =

1.0000 10.0000 24.0000 0

$\frac{Y(s)}{U(s)}=\frac{10s+20}{s^3+10s^2+24s}$

P=[-1+2j -1-2j -5]
K=acker(A,B,P)

>>

K =

1.2500 1.2500 0.1938

Aab=[1 0]
Abb=[0 1;-24 -10]
Pe=[-10 -10]
Ke=acker(Abb',Aab',Pe)'

>>

Ke =

10
-24

Aaa=[0]
Aba=[0;0]
Bb=[10;-80]
Ba=[0]
Ka=[1.25]
Kb=[1.25 0.1938]
Ahat=Abb-Ke*Aab
Bhat=Ahat*Ke+Aba-Ke*Aaa
Fhat=Bb-Ke*Ba
Atilde=Ahat-Fhat*Kb
Btilde=Bhat-Fhat*(Ka+Kb*Ke)
Ctilde=-Kb
Dtilde=-(Ka+Kb*Ke)
[num,den]=ss2tf(Atilde,Btilde,-Ctilde,-Dtilde)

>>

num =

9.0988 73.4880 125.0000

den =

1.0000 16.9960 -30.0400

$\frac{U(s)}{-Y(s)}=\frac{9.1s^2+73.488s+125}{s^2+16.996s-30.04}$

num=[10 20]
den=conv([1 0],conv([1 4],[1 6]))
sys1=tf(num,den)
num1=[9.1 73.5 125]
den1=[1 17 -30]
sys2=tf(num1,den1)
sys3=feedback(sys1,sys2,-1)

>>

sys3 =

10 s^3 + 190 s^2 + 40 s - 600
-------------------------------------------------
s^5 + 27 s^4 + 255 s^3 + 1025 s^2 + 2000 s + 2500

>> pzmap(sys3)

>> pzmap(sys2)

虽然闭环系统是稳定的，但由于内部控制器不稳定，所以此时选择的极点是不符合要求的，因此需要重新选择极点进行配置。

需要改变闭环极点的位置，或者观测器极点的位置，或者两者都要改变

Aab=[1 0]
Abb=[0 1;-24 -10]
Pe=[-4.5 -4.5]
Ke=acker(Abb',Aab',Pe)'

>>

Ke =

-1.0000
6.2500

Aaa=[0]
Aba=[0;0]
Bb=[10;-80]
Ba=[0]
Ka=[1.25]
Kb=[1.25 0.1938]
Ahat=Abb-Ke*Aab
Bhat=Ahat*Ke+Aba-Ke*Aaa
Fhat=Bb-Ke*Ba
Atilde=Ahat-Fhat*Kb
Btilde=Bhat-Fhat*(Ka+Kb*Ke)
Ctilde=-Kb
Dtilde=-(Ka+Kb*Ke)
[num,den]=ss2tf(Atilde,Btilde,-Ctilde,-Dtilde)

>>

num =

1.2112 11.2137 25.3125

den =

1.0000 5.9960 2.1295

$\frac{U(s)}{-Y(s)}=\frac{1.2s^2+11.2s+25}{s^2+6s+2.12}$

num=[10 20]
den=conv([1 0],conv([1 4],[1 6]))
sys1=tf(num,den)
num1=[1.2 11.2 25]
den1=[1 6 2.12]
sys2=tf(num1,den1)
pzmap(sys2)

sys3=feedback(sys1,sys2,-1)
pzmap(sys3)

此时系统不仅外部稳定而且内部也稳定

AA=[A-B*K B*Kb;zeros(2,3) Abb-Ke*Aab]
x0=[1 0 0 1 0]
t=0:0.01:10
x=initial(AA,eye(5),eye(5),eye(5),x0,t)
subplot(2,3,1)
plot(t,x(:,1)),grid,xlabel('t'),ylabel('x1')
subplot(2,3,2)
plot(t,x(:,2)),grid,xlabel('t'),ylabel('x2')
subplot(2,3,3)
plot(t,x(:,3)),grid,xlabel('t'),ylabel('x3')
subplot(2,3,4)
plot(t,x(:,4)),grid,xlabel('t'),ylabel('e1')
subplot(2,3,5)
plot(t,x(:,5)),grid,xlabel('t'),ylabel('e2')

比较：

1、状态反馈增益矩阵：

A=[0 1;0 -2]
B=[0;4]
C=[1 0]
D=0
P=[-2+2*sqrt(3)*j -2-2*sqrt(3)*j]
K=acker(A,B,P)

>>

K =

4.0000 0.5000

2、全维观测器增益矩阵

Pe=[-8 -8]
Ke=acker(A',C',Pe)'

>>

Ke =

14
36

零状态响应：

AA=[A-B*K B*K;zeros(2) A-Ke*C]
x0=[1 0 1 0]
t=0:0.01:10
x=initial(AA,eye(4),eye(4),eye(4),x0,t)
subplot(2,2,1)
plot(t,x(:,1)),grid,xlabel('t'),ylabel('x1')
subplot(2,2,2)
plot(t,x(:,2)),grid,xlabel('t'),ylabel('x2')
subplot(2,2,3)
plot(t,x(:,3)),grid,xlabel('t'),ylabel('e1')
subplot(2,2,4)
plot(t,x(:,4)),grid,xlabel('t'),ylabel('e2')

3、降维观测器增益矩阵

Aaa=[0]
Aab=[1]
Aba=[0]
Abb=[-2]
Ba=[0]
Bb=[4]
Pe=[-8]
Ke=acker(Abb',Aab',Pe)'

>>

Ke =

6

AAA=[A-B*K B*Kb;zeros(1,2) Abb-Ke*Aab]
t=0:0.01:10
x0=[1 0 1]
x=initial(AAA,eye(3),eye(3),eye(3),x0,t)
subplot(1,3,1)
plot(t,x(:,1)),grid,xlabel('t'),ylabel('x1')
subplot(1,3,2)
plot(t,x(:,2)),grid,xlabel('t'),ylabel('x2')
subplot(1,3,3)
plot(t,x(:,3)),grid,xlabel('t'),ylabel('e')

二者系统的bode图：

总结：

some optimization problems solved with matlab

t=0:0.01:8
for k=3:0.2:5
for a=0.1:0.1:3
num=[4*k 8*k*a 4*k*a^2]
den=[1 6 8+4*k 4+8*k*a 4*k*a^2]
y=step(num,den,t)
m=max(y)
if m<1.15&m>1.10
break 结束内循环
end
if m<1.15&m>1.10
break 结束外循环
end
end
end
solution=[k a m]
plot(t,y)
grid

>>

solution =

5.0000 0.7000 1.1019

?????????为什么我仿真出来的值和课本里面仿真出来的不一样

应该是m符合这个区间要求的k和a有很多组值。

k=3,a=1;m=1.1469

t=0:0.01:8
for k=[2:0.2:5]
for a=[0.1:0.1:3]
num=[4*k 8*k*a 4*k*a^2]
den=[1 6 8+4*k 4+8*k*a 4*k*a^2]
y=step(num,den,t)
m=max(y)
if m<2.6 & m>2.55
break
end
if m<2.6 & m>2.55
break
end
end
end
solution=[k a m]
plot(t,y)
grid

>>

solution =

5.0000 3.0000 2.9784

此时求出来的m不满足条件，因此计算出来的结果并不是真实的解

这是因为当条件不满足时，直到循环结束了才跳出来，因此取得最终值就是k和a的临界值。

for k=[2:0.2:3 NaN]
for a=[0.1:0.1:2 NaN]

>>

solution =

NaN NaN NaN

NaN:表示在这个区间内，没有值满足要求

因此在某个给定区间内，需要获得满足这一要求的所有的参数组的可行值。

t=0:0.01:8
for k=[3:0.2:5]
for a=[0.1:0.1:3]
num=[4*k 8*k*a 4*k*a^2]
den=[1 6 8+4*k 4+8*k*a 4*k*a^2]
y=step(num,den,t)
m=max(y)
if m<1.2 & m>1.1
break
end
if m<1.2 & m>1.1
break
end
end
end
solution=[k a m]
plot(t,y)
grid

>>

solution =

5.0000 0.7000 1.1019 这有一个解

K=3:0.2:5
a=0.1:0.1:3
k=0
for i=1:11
for j=1:30
num=[4*K(i) 8*K(i)*a(j) 4*K(i)*a(j)^2]
den=[1 6 8+4*K(i) 4+8*K(i)*a(j) 4*K(i)*a(j)^2]
y=step(num,den,t)
m=max(y)
if m<1.2&m>1.1
k=k+1
solution(k,:)=[K(i) a(j) m]
end
end
end

>> solution

solution =

3.0000 1.0000 1.1469
3.0000 1.1000 1.1979
3.2000 0.9000 1.1065
3.2000 1.0000 1.1581
3.4000 0.9000 1.1181
3.4000 1.0000 1.1688
3.6000 0.9000 1.1291
3.6000 1.0000 1.1791
3.8000 0.9000 1.1396
3.8000 1.0000 1.1888
4.0000 0.8000 1.1003
4.0000 0.9000 1.1497
4.0000 1.0000 1.1982
4.2000 0.8000 1.1107
4.2000 0.9000 1.1594
4.4000 0.8000 1.1208
4.4000 0.9000 1.1687
4.6000 0.8000 1.1304
4.6000 0.9000 1.1776
4.8000 0.8000 1.1396
4.8000 0.9000 1.1862
5.0000 0.7000 1.1019
5.0000 0.8000 1.1485
5.0000 0.9000 1.1945

sortsolution=sortrows(solution,3):对第三列，从小到大排列

t=0:0.01:8
k=2:0.1:3
a=0.5:0.1:1.5
k1=0
for i=1:11
for j=1:11
num1=conv([k(i)],conv([1 a(j)],[1 a(j)]))
den1=[1 0]
sys1=tf(num1,den1)
num2=[1.2]
den2=[0.36 1.86 2.5 1]
sys2=tf(num2,den2)
sys3=series(sys1,sys2)
sys4=feedback(sys3,1,-1)
y=step(sys4,t)
m=max(y)
if m<1.1
k1=k1+1
solution(k1,:)=[k(i) a(j) m]
end
end
end

>> solution

solution =

2.0000 0.5000 0.9552
2.0000 0.6000 0.9842
2.0000 0.7000 0.9958
2.0000 0.8000 0.9994
2.0000 0.9000 1.0614
2.1000 0.5000 0.9585
2.1000 0.6000 0.9856
2.1000 0.7000 0.9963
2.1000 0.8000 1.0077
2.1000 0.9000 1.0694
2.2000 0.5000 0.9614
2.2000 0.6000 0.9868
2.2000 0.7000 0.9966
2.2000 0.8000 1.0158
2.2000 0.9000 1.0772
2.3000 0.5000 0.9640
2.3000 0.6000 0.9878
2.3000 0.7000 0.9969
2.3000 0.8000 1.0239
2.3000 0.9000 1.0848
2.4000 0.5000 0.9663
2.4000 0.6000 0.9888
2.4000 0.7000 0.9972
2.4000 0.8000 1.0319
2.4000 0.9000 1.0923
2.5000 0.5000 0.9684
2.5000 0.6000 0.9896
2.5000 0.7000 0.9974
2.5000 0.8000 1.0398
2.5000 0.9000 1.0997
2.6000 0.5000 0.9703
2.6000 0.6000 0.9903
2.6000 0.7000 0.9976
2.6000 0.8000 1.0475
2.7000 0.5000 0.9719
2.7000 0.6000 0.9909
2.7000 0.7000 0.9978
2.7000 0.8000 1.0550
2.8000 0.5000 0.9735
2.8000 0.6000 0.9915
2.8000 0.7000 1.0024
2.8000 0.8000 1.0624
2.9000 0.5000 0.9749
2.9000 0.6000 0.9920
2.9000 0.7000 1.0101
2.9000 0.8000 1.0695
3.0000 0.5000 0.9761
3.0000 0.6000 0.9924
3.0000 0.7000 1.0177
3.0000 0.8000 1.0766

>> sortsolution=sortrows(solution,3)

sortsolution =

2.0000 0.5000 0.9552
2.1000 0.5000 0.9585
2.2000 0.5000 0.9614
2.3000 0.5000 0.9640
2.4000 0.5000 0.9663
2.5000 0.5000 0.9684
2.6000 0.5000 0.9703
2.7000 0.5000 0.9719
2.8000 0.5000 0.9735
2.9000 0.5000 0.9749
3.0000 0.5000 0.9761
2.0000 0.6000 0.9842
2.1000 0.6000 0.9856
2.2000 0.6000 0.9868
2.3000 0.6000 0.9878
2.4000 0.6000 0.9888
2.5000 0.6000 0.9896
2.6000 0.6000 0.9903
2.7000 0.6000 0.9909
2.8000 0.6000 0.9915
2.9000 0.6000 0.9920
3.0000 0.6000 0.9924
2.0000 0.7000 0.9958
2.1000 0.7000 0.9963
2.2000 0.7000 0.9966
2.3000 0.7000 0.9969
2.4000 0.7000 0.9972
2.5000 0.7000 0.9974
2.6000 0.7000 0.9976
2.7000 0.7000 0.9978
2.0000 0.8000 0.9994
2.8000 0.7000 1.0024
2.1000 0.8000 1.0077
2.9000 0.7000 1.0101
2.2000 0.8000 1.0158
3.0000 0.7000 1.0177
2.3000 0.8000 1.0239
2.4000 0.8000 1.0319
2.5000 0.8000 1.0398
2.6000 0.8000 1.0475
2.7000 0.8000 1.0550
2.0000 0.9000 1.0614
2.8000 0.8000 1.0624
2.1000 0.9000 1.0694
2.9000 0.8000 1.0695
3.0000 0.8000 1.0766
2.2000 0.9000 1.0772
2.3000 0.9000 1.0848
2.4000 0.9000 1.0923
2.5000 0.9000 1.0997

k=2
a=0.7
num1=conv([k],conv([1 a],[1 a]))
den1=[1 0]
sys1=tf(num1,den1)
num2=[1.2]
den2=[0.36 1.86 2.5 1]
sys2=tf(num2,den2)
sys3=series(sys1,sys2)
sys4=feedback(sys3,1,-1)
y=step(sys4,t)
plot(t,y)
grid

t=0:0.01:8
k=3:0.2:5
a=0.1:0.2:3
k1=0
for i=1:11
for j=1:15
num1=conv([k(i)],conv([1 a(j)],[1 a(j)]))
den1=[1 0]
sys1=tf(num1,den1)
num2=[1.2]
den2=[0.36 1.86 2.5 1]
sys2=tf(num2,den2)
sys3=series(sys1,sys2)
sys4=feedback(sys3,1,-1)
y=step(sys4,t)
m=max(y)
s=801
while y(s)>0.98&y(s)<1.02 调节时间是怎么定义的
s=s-1
end
ts=(s-1)*0.01 y(s):引用？？？？为什么要乘以0.01
if m<1.15&m>1.10
k1=k1+1
solution(k1,:)=[k(i) a(j) m ts]
end
end
end

k=4.4
a=0.7
num1=conv([k],conv([1 a],[1 a]))
den1=[1 0]
sys1=tf(num1,den1)
num2=[1.2]
den2=[0.36 1.86 2.5 1]
sys2=tf(num2,den2)
sys3=series(sys1,sys2)
sys4=feedback(sys3,1,-1)
y=step(sys4,t)
plot(t,y)
grid

Quadratic optimal control：二次型的最优控制

问题：

A=[-1 1;0 2]
B=[1;0]
Q=[1 0;0 1]
R=[1]
K=lqr(A,B,Q,R)

>>

错误使用 lqr (第 42 行)
Cannot compute a stabilizing LQR gain (the Riccati solution S and gain matrix K are
infinite).
This could be because:
* A has unstable modes that are not controllable through B,
* Q,R,N values are too large,
* [Q N;N' R] is indefinite,
* The E matrix in the state equation is singular.

出错 main (第 7 行)
K=lqr(A,B,Q,R)

syms k1 k2 s
A=[-1 1;0 2]
B=[1;0]
k=[k1 k2]
I=eye(2)
det(s*I-A+B*k)

>>

ans =

(s - 2)*(k1 + s + 1)

不论k怎么去，闭环你系统始终有一个极点在s的右半平面，即系统始终不稳定。

因此这个系统不能使用最优控制

A=[0 1 0;0 0 1;-35 -27 -9]
B=[0;0;1]
Q=[1 0 0;0 1 0;0 0 1]
R=[1]
[K,P,E]=lqr(A,B,Q,R)

>>

K =

0.0143 0.1107 0.0676

P =

4.2625 2.4957 0.0143
2.4957 2.8150 0.1107
0.0143 0.1107 0.0676

E =

-1.9859 + 1.7110i
-1.9859 - 1.7110i
-5.0958 + 0.0000i

AA=A-B*K
x0=[1;0;0]
t=0:0.1:10
eye(3,1)
[y,x]=initial(AA,eye(3,1),eye(1,3),0,x0,t)
subplot(1,3,1)
plot(t,x(:,1))
subplot(1,3,2)
plot(t,x(:,2))
subplot(1,3,3)
plot(t,x(:,3))

最优控制就相当于对于系统性能的评价方面，选定一个标量函数，当该标量函数最小的时候，就可以实现最优控制了，但实现的方法是通过状态反馈来实现的。

A=[0 1 0;0 0 1;0 -2 -3]
B=[0;0;1]
C=[1 0 0]
D=[0]
Q=[100 0 0;0 1 0;0 0 1]
R=[0.01]
[k,P,E]=lqr(A,B,Q,R)

>>

k =

100.0000 53.1200 11.6711

P =

55.1200 14.6711 1.0000
14.6711 7.0267 0.5312
1.0000 0.5312 0.1167

E =

-2.2141 + 2.2047i
-2.2141 - 2.2047i
-10.2429 + 0.0000i

AA=A-B*k
BB=100*B
t=0:0.01:10
[y,x,t]=step(AA,BB,C,D,1,t)
subplot(2,2,1)
plot(t,x(:,1))
subplot(2,2,2)
plot(t,x(:,2))
subplot(2,2,3)
plot(t,x(:,3))
subplot(2,2,4)
plot(t,y)

Appendix

num=[1 3]
den=[1 5 20 16 0]
rlocus(num,den)
v=[-6 6 -6 6]
axis(v)
grid

grid

可以用两次grid命令去除栅格线

x1=[-6 4]
y1=[0 0]
plot(x1,y1)

hold on
x2=[0 0]
y2=[-6 6]
plot(x2,y2)

x=[1 1 1 1 1 4 4 4 4 4 ]
y=[2 1 0 -1 -2 -2 -1 0 1 2]
plot(x,y)
v=[0 5 -3 3]
axis(v)