| |
lsimは、コントロールツールボックスの中の関数のためコントロールツールボックスを購入しないと動作しない．
しかし，ない場合には，自分で互換の関数を作ってしまうことも可能である．ここでは，
参考としてlsim互換のスクリプトプログラムを挙げておく．
7通りの伝達関数サンプルでテストし，ほぼlsimと同じ結果が出ているので問題はないと思います。
-------------lsim.m互換プログラム-------llsim.m---
% ご指摘で解ったのですが、lsimやc2dmはコントロールツールボックスに付属する
% Mファイルのため購入する必要があります。
%このプログラムは、MATLAB Handbook のP.100で使用しているlsim互換プログラム
% です。
% MATLABのみの標準的な関数で伝達関数のシミュレーションします。
% 解き方は、まず伝達関数を状態方程式行列に分解します。
% 出来た状態方程式をPade'近似により離散系に変換し解いています。
%
%
% 入力形式は、[Y,T]=llsim([分子],[分母],[入力信号],[時間]);
% となっています。
%
% Sample Program
%
% tau = 0.01;
%[Y,TT]=llsim(1,[1 1],ones(size(0:tau:10)),0:tau:10);
%[Y,TT]=llsim(1,conv([T 1],[T 1]),ones(size(0:tau:10)),0:tau:10);
%のように入力します。
% Programmed by Kazuyuki Kobayashi, 1998.5.18
% E-mail: ikko@k.hosei.ac.jp
% 動作確認は、MATLAB Version 4.2c,5.2で行ないました。
%
function [Y,T]=llsim(NUM,DEN,U,T)
%参考文献: システム制御:増淵正美,コロナ社，pp.70-71
A=zeros(length(DEN)-1);
DEN1 = DEN(1);
DEN = DEN / DEN1;
NUM = NUM / DEN1;
% NUMとDENのサイズを同じにする。
shiftnum = length(DEN)-length(NUM);
NUM = [zeros(1,shiftnum) NUM];
%式変形したベータ（分子）を計算する
BNUM=zeros(size(DEN));
for ii=1:length(DEN)
BNUM(ii)=NUM(ii)-NUM(1)*DEN(ii);
end;
%伝達関数を可制御標準形の状態方程式に直す。
%
% define A matrix
if length(DEN) > 2
for ii=1:length(DEN)-2
A(ii,1+ii)=1;
end
A(length(DEN)-1,:)=-fliplr(DEN(2:length(DEN)));
else
A=-DEN(2);
end
% define B matrix
if length(DEN) > 1
B=zeros(length(DEN)-1,1);
B(length(DEN)-1,1)=1;
else B = 1;
end
%define C and D matrix
if(length(BNUM) > 1)
D=NUM(length(NUM(1)));
C=fliplr(BNUM(2:length(BNUM)));
else
C=BBNUM(length(BBNUM));
D=0;
end
% 連続系の状態方程式をPade'近似に離散系に変換する。
% 参考文献:制御工学,渡辺嘉二郎著，サイエンス社
tau = T(2)-T(1);
[raw col]= size(A);
n= col; % X
[raw col]= size(C);
m= col; % U
I=eye(n);
H=tau * inv(I - 0.5 * tau * A);
P = A * H + I;
Q = H * B;
Y=zeros(size(T));
X=zeros(n,1);
for ii=1:length(T)
X = P * X + Q * U(1,ii);
Y(ii) = C * X + D * U(1,ii);
end
------------以降は、検証用に使用したプログラムです。----------------
clear
close all
T=1;
tau = 0.01;
TTT=0:tau:10;
Y=llsim(1,[T 1],ones(size(TTT)),TTT);
YY=lsim(1,[T 1],ones(size(TTT)),TTT);
figure;subplot(2,1,1);plot(TTT,Y,TTT,YY);
subplot(2,1,2);plot(TTT,Y)
Y=llsim(1,conv([T 1],[T 1]),ones(size(TTT)),TTT);
YY=lsim(1,conv([T 1],[T 1]),ones(size(TTT)),TTT);
figure;subplot(2,1,1);plot(TTT,Y,TTT,YY);
subplot(2,1,2);plot(TTT,Y);
Y=llsim(1,conv([T 1],conv([T 1],[T 1])),ones(size(TTT)),TTT);
YY=lsim(1,conv([T 1],conv([T 1],[T 1])),ones(size(TTT)),TTT);
figure;subplot(2,1,1);plot(TTT,Y,TTT,YY);
subplot(2,1,2);plot(TTT,Y)
Y=llsim([T 0],[T 1],ones(size(TTT)),TTT);
YY=lsim([T 0],[T 1],ones(size(TTT)),TTT);
figure;subplot(2,1,1);plot(TTT,Y,TTT,YY);
subplot(2,1,2);plot(TTT,Y)
YY=lsim([1 1],[1 2 1],ones(size(TTT)),TTT);
Y=llsim([1 1],[1 2 1],ones(size(TTT)),TTT);
figure;subplot(2,1,1);plot(TTT,Y,TTT,YY);
subplot(2,1,2);plot(TTT,Y)
YY=lsim([1 2],[1 2 1],ones(size(TTT)),TTT);
Y=llsim([1 2],[1 2 1],ones(size(TTT)),TTT);
figure;subplot(2,1,1);plot(TTT,Y,TTT,YY);
subplot(2,1,2);plot(TTT,Y)
YY=lsim([1 2 1],[1 2 2 1],ones(size(TTT)),TTT);
Y=llsim([1 2 1],[1 2 2 1],ones(size(TTT)),TTT);
figure;subplot(2,1,1);plot(TTT,Y,TTT,YY);
subplot(2,1,2);plot(TTT,Y)
end
|