// pid_lagdelay.sce		2009/08/01
// 1次遅れ＋むだ時間系（目標値応答）
s=%s;
// プラントのパラメータ
K=0.98;							// ゲイン
T=4.85;							// 時定数[s]
L=4;							// むだ時間[s]
G=K/(T*s+1)*(12-6*L*s+L^2*s^2)/(12+6*L*s+L^2*s^2);	// プラントの伝達関数（むだ時間を2次パデ近似）
G1=K*(12-6*L*s+L^2*s^2)/(12+6*L*s+L^2*s^2);		// むだ時間系
G2=G1/s;						// PI補償器設計用の仮プラント

// PI補償器の設計
Kpt=input('Enter Kpt');				// 設計パラメータKp_tildeの入力（Kpt-->大 のとき速応性-->大）
L=Kpt*G2;
[p,fr]=p_margin(L);
wp=fr*2*%pi// ゲイン交差周波数
pm=180+p// 位相余裕
phi_M=input('Enter phi_M');			// 設計パラメータの位相余裕（°）の入力
phia=(phi_M-pm)*%pi/180;			// 必要な位相進み量
if phia > 0 then				// 必要な位相進み量を与えるT_Iの計算
  T_I=tan(phia)/wp;
else
  T_I=0;
end
Kp1=Kpt/abs(1+%i*wp*T_I)			// ゲインの補正

// カスケード形式のPIDパラメータを標準的なPIDパラメータに置き換える。
kappa_P=Kp1*(T_I+T);
tau_I=T_I+T;
tau_D=T_I*T/(T_I+T);

eps=0.1;
s1=s/(1+eps*tau_D*s);				// 近似微分器
C=kappa_P*(1+1/(tau_I*s)+tau_D*s1);		// PID補償器
C1=kappa_P*(1+1/(tau_I*s));			// 偏差eに対する比例＋積分動作
C2=kappa_P*tau_D*s1;				// 出力フィードバックに対する微分動作

Lo=G*C;// ループ伝達関数
[p,fr]=p_margin(Lo);
wp=fr*2*%pi					// ゲイン交差周波数
pm=180+p					// 位相余裕

// 閉ループ伝達関数の計算
L1=C1*G/(1+G*C2);
Gyr=L1/(1+L1);					// 目標入力rから出力yまでの伝達関数
Gyd=G/(1+Lo);					// 外乱dからyまでの伝達関数
Gur=C1/(1+Lo);					// rからuまでの伝達関数
Gud=-Lo/(1+Lo);					// dからuまでの伝達関数

// 連続時間線形系として設定
syr=syslin('c',Gyr);
syd=syslin('c',Gyd);
sur=syslin('c',Gur);
sud=syslin('c',Gud);

// ステップ応答のシミュレーション
t=0:0.01:40;
yr=csim('step',t,syr);
yd=csim('step',t,syd);
ur=csim('step',t,sur);
ud=csim('step',t,sud);

// 応答のプロット
xset('window',0)
plot2d(t,ur,6)
xgrid(4)

xset('window',1)
plot2d(t,yr,6)
xgrid(4)

xset('window',4)
plot2d(t,ud,6)

xset('window',2)
xgrid(4)
plot2d(t,yd,6)

// 設計されたPIDパラメータ
kappa_P
tau_I
tau_D

// ゲイン交差周波数と位相余裕
wp
pm