MATLAB离散最优控制程序应该怎么编？

1、step的调用方式有问题，完全没有体现出这是一个离散系统。考虑：

step(ss(Ac,Bc,Cc,Dc,1))

2、【k1=k(1)Ac=A-B*kBc=B*k1Cc=CDc=D】Bc和Cc的计算有问题，应该是Bc=B, Cc=C-D*k。

3、LQR最优调节器与Q、R的选择关系很大，所谓的最优只是对于特定的Q、R而言的。

4、可以试一试用下面的代码观察各状态量：

step(ss(Ac,Bc,eye(length(A)),0,1))

一、连续Lyapunov方程

连续Lyapunov方程可以表示为

Lyapunov方程来源与微分方程稳定性理论，其中要求C为对称正定的n×n方阵，从而可以证明解X亦为n×n对称矩阵，这类方程直接求解比较困难，不过有了Matlab中lyap()函数，就简单多了。

>>A=[1 2 34 5 67 8 0]

A =

1 2 3

4 5 6

7 8 0

>>C=-[10 5 45 6 74 7 9]

C =

-10-5-4

-5-6-7

-4-7-9

>>X=lyap(A,C)

X =

-3.94443.88890.3889

3.8889 -2.77780.2222

0.38890.2222 -0.1111

二、Lyapunov方程的解析解

利用Kroncecker乘积的表示方法，可以将Lyapunov方程写为

function x=lyap2(A,C)

%Lyapunov方程的符号解法

n=size(C,1)

A0=kron(A,eye(n))+kron(eye(n),A)

c=C(:)

x0=-inv(A0)*c

x=reshape(x0,n,n)

例子

>>A=[1 2 34 5 67 8 0]

>>C=-[10 5 45 6 74 7 9]

>>x=lyap2(sym(A),sym(C))

x =

[ -71/18, 35/9, 7/18]

[ 35/9, -25/9,2/9]

[ 7/18,2/9, -1/9]

三、离散Lyapunov方程

离散Lyapunov方程的一般形式为

Matlab中直接提供了dlyap()函数求解该方程：X=dlyap(A,Q)

其实，如果A矩阵非奇异，则等式两边同时右乘得到

就可以将其变换成连续的Sylvester方程

而Sylvester方程是广义Lyapunov方程，故离散的Lyapunov方程还可以使用下面的方法求解

B=-inv(A’)

C=Q*inv(A’)

X=lyap(A,B,C)

下面总结下我们上面的讲到的知识点：

X=lyap(A,C) 连续Lyapunov方程数值解法

X=lyap2(A,C) 连续Lyapunov方程符号解法

X=lyap(A,B,C)广义Lyapunov方程，即Sylvester方程

X=dlyap(A,Q)或者X=lyap(A,-inv(A’),Q*inv(A’))离散Lyapunov方程

Sylvester方程Matlab求解Sylvester方程的一般形式为

该方程又称为广义的Lyapunov方程，式中A是n×n方阵，B是m×m方阵，X和C是n×m矩阵。Matlab控制工具箱提供了直接的求解该方程的lyap()函数

A=[8 1 63 5 74 9 2]

B=[2 34 5]

C=[1 23 45 6]

X=lyap(A,B,C)

A =

8 1 6

3 5 7

4 9 2

B =

2 3

4 5

C =

1 2

3 4

5 6

X =

0.20110.2016

0.03930.1554

-0.6428 -0.8966

同理，我们使用Kronecker乘机的形式将原方程进行如下变化

故可以编写Sylvester方程的解析解函数

function X=lyap3(A,B,C)

%Sylvester方程的解析解法

%rewrited by dynamic

%more information http://www.ilovematlab.cn

If nargin==2,C=BB=A'end

[nr,nc]=size(C)

A0=kron(A,eye(nc))+kron(eye(nr),B')

try

C1=C'

X0=-inv(A0)*C1(:)

X=reshape(X0,nc,nr)

catch

error('Matlabsky提醒您：矩阵奇异！')

end

用上面的数据，我们试验下该解析解法的

>>X=lyap3(sym(A),B,C)

X =

[ 2853/14186,557/14186, -9119/14186]

[ 11441/56744, 8817/56744, -50879/56744]

Riccati方程的Matlab求解Riccati方程是一类很著名的二次型矩阵形式，其一般形式为

由于含有矩阵X的二次项，所有Riccati方程求解要Lyapunov方程更难，Matlab控制工具箱提供了are()函数，可以直接求解该函数

A=[-2 1 -3-1 0 -20 -1 -2]

B=[2 2 -2-1 5 -2-1 1 2]

C=[5 -4 41 0 41 -1 5]

X=are(A,B,C)

A =

-2 1-3

-1 0-2

0-1-2

B =

2 2-2

-1 5-2

-1 1 2

C =

5-4 4

1 0 4

1-1 5

X =

0.9874 -0.79830.4189

0.5774 -0.13080.5775

-0.2840 -0.07300.6924

如何用matlab求解lyapunov指数我是需要分析计算LOGISTIC数据 ,都是用来说明对初值的敏感.以下是LOGISTIC求解的程序 ,希望得到LYAPUNOV的程序.

clc

clear

close all

lambda = 3:5e-4:4

x = 0.4*ones(1,length(lambda))

N1 = 400 % 前面的迭代点数

N2 = 100 % 后面的迭代点数

f = zeros(N1+N2,length(lambda))

for i = 1:N1+N2

x = lambda .* x .* (1 - x)

f(i,:) = x

end

f = f(N1+1:end,:)

plot(lambda,f,'k.','MarkerSize',1)

xlabel('\lambda')

ylabel('x')

慢慢看吧，很有用

---