ここでは、次の４種類の積分動作をもつスライディングモード制御（SMI制御）系の設計について述べます。

Note e13に示すように、データセットNo.1について、id=128の場合を設計用モデルとして採用したとき、Case2とCase3のどちらを用いても制御性能に大きな差異が見られません。そこでこのコントローラがどれくらいのモデル変動をカバーできるかという意味のロバスト性をシミュレーションにより調べてみます。ここで、モデル変動としては、前進速度を変える場合と、重心位置を変える場合の２通りを考えます。そのために次の３組のデータセットを考えます。

設計用モデルとして、データセットNo.3における船外機取付角-2[deg]のときのid=128の場合を採用し、Case 2とCase 3のLQIコントローラを設計します。そしてデータセットNo.3,4,5における船外機取付角-2[deg]のときのモデル（id=126～133）を安定化できるかを調べます。シミュレーション結果を次に示します。

⇒Case2:いくつかの動作点で不安定
⇒Case3:すべて安定
図1　データセットNo.3の場合、SMIコントローラのロバスト性

⇒Case2:いくつかの動作点で不安定
⇒Case3:すべて安定
図2　データセットNo.4の場合、SMIコントローラのロバスト性

⇒Case2:いくつかの動作点で不安定
⇒Case3:すべて安定
図3　データセットNo.5の場合、SMIコントローラのロバスト性

これらのシミュレーションは本来は非線形シミュレータ上で行うべきですが、ここでは６次元線形モデルを用いて簡易的に行っています。

第１の注意点は、６次元線形モデルの第４番目の状態変数の係数を強制的にとしていることです。これはとなる場合があって、が発散してしまうからです。

第２の注意点は、上のシミュレーションではCase 2のコントローラではすべてを安定化できていませんが、設計用モデルとして、データセットNo.1における船外機取付角-2[deg]のときのid=128の場合を採用すると、すべて安定化されることです。

これらの事情は、線形モデルが暫定版の非線形シミュレータから得られているので正確ではないためとも考えられます。いずれにしてもできるだけ実機のダイナミックスを反映した非線形シミュレータ上での再検討が必要になります。

%cYMCPB5_robust_smi.m
%-----
 clear all, close all
 dataset=[...
 "ABC_No1_20220128",...
 "ABC_No2_20220131",...
 "ABC_No3_20220128",...
 "ABC_No4_20220131",...
 "ABC_No5_20220131"];
%-----設計用データセット
 no=3
 load(dataset(no))
 ID=[7,8,9,10,11,12,13,14,...
 36,37,38,39,40,41,42,43,...
 66,67,68,69,70,71,72,73,...
 96,97,98,99,100,101,102,103,...
 126,127,128,129,130,131,132,133];
%-----６次元モデル
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 sys=ss(A0(:,:,ID),B0(:,:,ID),eye(6,6),[]);
%-----５次元モデル
%(1)z,(2)th,(3)dz,(4)dth,(5)the
 k=[2,3,5,6]; i=[2,4,5]; j=[2];
 for id=ID
   A(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,5)];
 end 
 omegae=0.0380*3;
 B=[zeros(4,1);omegae];
 E=eye(5,5);
 CM=E(i,:);
 C=CM(3,:);
 sys5=ss(A(:,:,ID),B,eye(5,5),[]);  
%-----設計ポイント 
 index=35 %[3,11,19,27,35]
 id=ID(index); 
 Vs=dx_eq_Series(id);    %*3.6
 zs=z_eq_Series(id);     %*100
 ths=th_eq_Series(id);   %/pi*180
 thes=the_eq_Series(id); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,index));
 [A,B]=ssdata(sys5(:,:,index));
%----- 
% Commands to make a unit vec tor controller incorporating integral action.
% It assumes the triple (A,B,C) is already in regular form
% For simplicity robust pole placement is used for the hyperplane design
% (with suitable modification any of the methods in chap 4 could be used)
% the nth order vector p1 specifies the reduced order motion poles
% the mth order vector p2 specifies  robust pole placement
% the range space dynamics
%-----
 [nn,mm]=size(B);
 [pp,nn]=size(C);
% Augment the orignal (A,B) with the integrator states
 bigA=[zeros(pp,pp) -C;
       zeros(nn,pp) A];
 bigB=[zeros(pp,mm);
       B];
% Since it is assumed that (A,B) is in regular form the augmented 
% pair (bigA,bigB) is also automatically in refular form
 A11=bigA(1:pp+nn-mm,1:pp+nn-mm);
 A12=bigA(1:pp+nn-mm,nn+pp-mm+1:nn+pp);
 B2=bigB(nn+pp-mm+1:nn+pp,:);
% Hyperplane design using robust pole placement
 j=sqrt(-1);
 p1=[-4.7/20+j*4.7, -4.7/20-j*4.7, -7.8/20+j*7.8, -7.8/20-j*7.8, -2/10];
 M=place(A11,A12,p1);
 A11s=A11-A12*M;
 S2=eye(mm)  % For simplicity
 S=S2*[M eye(mm)];
% Form the range space dynamics system matrix
 p2=-1
 Phi=diag(p2,0);
 P2hat=lyap(Phi,eye(mm));
 P2=inv(P2hat);
% Calculate the time varying component of the switching function
 Br=[eye(pp);
     zeros(nn-mm,pp)];
 Sr=-S2*inv(Br'*inv(A11s)*A12)*Br'*inv(A11s)*Br;
% Computes the gains for the feedback and feed-forward components
 L=-inv(S*bigB)*S*bigA+inv(S*bigB)*Phi*S
 L(2)=0; L(4)=0;
 Lr=-inv(S*bigB)*(Phi*Sr+S(:,1:pp))
 Ldr=inv(S*bigB)*Sr
 rho=1
 Ln=inv(S*bigB)*rho
%-----全モデルに対するシミュレーション
for no=[3,4,5] 
for iSIM=index-2:index+5
 id=ID(iSIM);
 Vs=dx_eq_Series(id);    %*3.6
 zs=z_eq_Series(id);     %*100
 ths=th_eq_Series(id);   %/pi*180
 thes=the_eq_Series(id); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,iSIM));
 Asys(4,4)=0;
 [A,B]=ssdata(sys5(:,:,iSIM));  
 us=thes;
 xs=[0;zs;ths;Vs;0;0];
 x0=[0;0;-1/180*pi;0;0;0];
 yc=0; 
 E=eye(6,6);
 CM=E([2,3,5,6],:);
 CM1=100*E([2],:);
 CM2=180/pi*E([3],:);
 ns=4;
 simout=sim('nYMCPB_smi.slx');
%[no,iSIM,id,180/pi*us]
 figure(no)
 subplot(221),plot(simout.x.Time,100*simout.x.Data(:,2)),hold on
 subplot(223),plot(simout.y.Time,180/pi*simout.y.Data(:,1)),hold on
 subplot(222),plot(simout.y.Time,180/pi*simout.y.Data(:,3)),hold on
 subplot(224),plot(simout.u.Time,simout.u.Data),hold on
end %iSIM
 figure(no)
 subplot(221),title("z[cm]"),grid
 subplot(223),title("th[deg]"),grid
 subplot(222),title("the[deg]"),grid
 subplot(224),title("ue"),grid
end %no 
%-----
%eof

%cYMCPB3_robust_smi.m
%-----
 clear all, close all
 dataset=[...
 "ABC_No1_20220128",...
 "ABC_No2_20220131",...
 "ABC_No3_20220128",...
 "ABC_No4_20220131",...
 "ABC_No5_20220131"];
%-----設計用データセット
 no=3
 load(dataset(no))
 ID=[7,8,9,10,11,12,13,14,...
 36,37,38,39,40,41,42,43,...
 66,67,68,69,70,71,72,73,...
 96,97,98,99,100,101,102,103,...
 126,127,128,129,130,131,132,133];
%-----６次元モデル
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 sys=ss(A0(:,:,ID),B0(:,:,ID),eye(6,6),[]);
%-----３次元モデル
%(1)th,(2)dth,(3)the
 k=[3,6]; j=[2];
 for id=ID
   AA(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,3)];
 end 
 omegae=0.0380*3;
 B=[zeros(2,1);omegae];
 E=eye(3,3);
 C=E(3,:);
 sys3=ss(AA(:,:,ID),B,eye(3,3),[]);  
%-----設計ポイント 
 index=35 %[3,11,19,27,35]
 id=ID(index); 
 Vs=dx_eq_Series(id);    %*3.6
 zs=z_eq_Series(id);     %*100
 ths=th_eq_Series(id);   %/pi*180
 thes=the_eq_Series(id); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,index));
 [A,B]=ssdata(sys3(:,:,index));
%----- 
% Commands to make a unit vec tor controller incorporating integral action.
% It assumes the triple (A,B,C) is already in regular form
% For simplicity robust pole placement is used for the hyperplane design
% (with suitable modification any of the methods in chap 4 could be used)
% the nth order vector p1 specifies the reduced order motion poles
% the mth order vector p2 specifies  robust pole placement
% the range space dynamics
%-----
 [nn,mm]=size(B);
 [pp,nn]=size(C);
% Augment the orignal (A,B) with the integrator states
 bigA=[zeros(pp,pp) -C;
       zeros(nn,pp) A];
 bigB=[zeros(pp,mm);
       B];
% Since it is assumed that (A,B) is in regular form the augmented 
% pair (bigA,bigB) is also automatically in refular form
 A11=bigA(1:pp+nn-mm,1:pp+nn-mm);
 A12=bigA(1:pp+nn-mm,nn+pp-mm+1:nn+pp);
 B2=bigB(nn+pp-mm+1:nn+pp,:);
% Hyperplane design using robust pole placement
 j=sqrt(-1);
 p1=[-4.7/20+j*4.7, -4.7/20-j*4.7, -2/10];
 M=place(A11,A12,p1);
 A11s=A11-A12*M;
 S2=eye(mm)  % For simplicity
 S=S2*[M eye(mm)];
% Form the range space dynamics system matrix
 p2=-1
 Phi=diag(p2,0);
 P2hat=lyap(Phi,eye(mm));
 P2=inv(P2hat);
% Calculate the time varying component of the switching function
 Br=[eye(pp);
     zeros(nn-mm,pp)];
 Sr=-S2*inv(Br'*inv(A11s)*A12)*Br'*inv(A11s)*Br;
% Computes the gains for the feedback and feed-forward components
 L=-inv(S*bigB)*S*bigA+inv(S*bigB)*Phi*S
 Lr=-inv(S*bigB)*(Phi*Sr+S(:,1:pp))
 Ldr=inv(S*bigB)*Sr
 rho=1
 Ln=inv(S*bigB)*rho 
%-----全モデルに対するシミュレーション
for no=[3,4,5] 
for iSIM=index-2:index+5
 id=ID(iSIM);
 Vs=dx_eq_Series(id);    %*3.6
 zs=z_eq_Series(id);     %*100
 ths=th_eq_Series(id);   %/pi*180
 thes=the_eq_Series(id); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,iSIM));
 Asys(4,4)=0;
 [A,B]=ssdata(sys3(:,:,iSIM));    
 us=thes;
 xs=[0;zs;ths;Vs;0;0];
 x0=[0;0;-1/180*pi;0;0;0];
 yc=0; 
 E=eye(6,6);
 CM=E([3,6],:);
 CM1=100*E([2],:);
 CM2=180/pi*E([3],:);
 ns=2;
 simout=sim('nYMCPB_smi.slx');
%[no,iSIM,id,180/pi*us]
 figure(no)
 subplot(221),plot(simout.x.Time,100*simout.x.Data(:,2)),hold on
 subplot(223),plot(simout.y.Time,180/pi*simout.y.Data(:,1)),hold on
 subplot(222),plot(simout.y.Time,180/pi*simout.y.Data(:,3)),hold on
 subplot(224),plot(simout.u.Time,simout.u.Data),hold on
end %iSIM
 figure(no)
 subplot(221),title("z[cm]"),grid
 subplot(223),title("th[deg]"),grid
 subplot(222),title("the[deg]"),grid
 subplot(224),title("ue"),grid
end %no  
%-----
%eof

Note e13　SMI制御系の設計
Case 1：５次元モデルに基づくSMI制御

●制御対象

の出力を、コマンド（次式の解）

に追従させることを考えます（は安定行列）。そのために、積分動作

を導入し、次の拡大系を構成します。ここで、(1)はすでに標準形であるとしています。

これを、次のように分割し直しても標準形であることには変わりありません。

ただし

●この積分器による拡大系を安定化できれば、積分器の値は定値となり、被積分項の値は零となり、はへ漸近します。そこで、SM制御によって拡大系を安定化し、追従制御系を構成することを考えます。この制御系は特別なの場合を含みますので、まずスイッチング関数として、次式を考えます。

(5)に対して、座標変換

を行って、次式を得ます。

ただし

ここで、が安定行列となるように行列が選ばれているとします。

●このときSM制御則は次式で与えられます。

%cYMCPB5_smi.m
%-----
 clear all, close all
 dataset=[...
 "ABC_No1_20220128",...
 "ABC_No2_20220131",...
 "ABC_No3_20220128",...
 "ABC_No4_20220131",...
 "ABC_No5_20220131"];
 no=1
 load(dataset(no))
 ID=[8,9,10,11,12,13,14,...
 36,37,38,39,40,41,42,43,...
 66,67,68,69,70,71,72,73,...
 97,98,99,100,101,102,103,...
 126,127,128,129,130,131,132,133];
%-----６次元モデル
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 sys=ss(A0(:,:,ID),B0(:,:,ID),eye(6,6),[]);
%-----５次元モデル
%(1)z,(2)th,(3)dz,(4)dth,(5)the
 k=[2,3,5,6]; i=[2,4,5]; j=[2];
 for id=ID
   A(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,5)];
 end 
 omegae=0.0380*3;
 B=[zeros(4,1);omegae];
 E=eye(5,5);
 CM=E(i,:);
 C=CM(3,:);
 sys5=ss(A(:,:,ID),B,eye(5,5),[]);  
%-----設計ポイント id=128=ID(33)
 Vs=dx_eq_Series(128);    %*3.6
 zs=z_eq_Series(128);     %*100
 ths=th_eq_Series(128);   %/pi*180
 thes=the_eq_Series(128); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,33));
 [A,B]=ssdata(sys5(:,:,33));
%-----
 AE=[0 -C ;
    zeros(5,1) A];
 BE=[0; B];
 lambda=[-4.7/20 4.7 -7.8/20 7.8 -2/10];
 nocomp=2;
 specpos=rand(6,5);
 specent=rand(6,5);
 S=swfvpl(AE,BE,lambda,nocomp,specpos,specent);   
 FE=(S*BE)\(S*AE);
 eig(AE-BE*FE);
%-----
 [nn,mm]=size(B);
 S1=S(:,1:nn-mm);
 S2=S(:,nn-mm+1:nn);
 Br=eye(nn-mm,mm);
 B2=B(nn-mm+1:nn,:); 
 Lambda=S2*B2;
 Phi=-0.5
 P2=0.5*inv(-Phi);
 Leq=inv(S*BE)*S*AE;
 LPhi=-inv(S*BE)*Phi*S;
 rho=2
 Sr=0
%----- 
 L=Leq+LPhi
 Ln=inv(S*BE)*rho
 Lr=inv(Lambda)*(Phi*Sr+S1*Br)
 Ldr=inv(Lambda)*Sr
%L(2)=0; L(4)=0;
%-----
 us=thes;
 xs=[0;zs;ths;Vs;0;0];
 x0=[0;0;-1/180*pi;0;0;0];
 yc=0; 
 E=eye(6,6);
 CM=E([2,3,5,6],:);
 CM1=100*E([2],:);
 CM2=180/pi*E([3],:);
 ns=4;
%-----
%eof
% L =
%    -0.6767   53.6065  -82.8461    3.7131   18.0214    2.8876
% Ln =
%    39.6861
% Lr =
%   -11.3251
% Ldr =
%      0

図1　SMI制御系シミュレーション

図2　５次元モデルに基づくSMI制御系と符号関数を通したSMI制御

Case 2：Case1でヒーブゲインを零とする場合

図3　ヒーブゲインを零とした場合のSMI制御系と符号関数を通したSMI制御

Case 3：３次元モデルに基づくSMI制御

●制御対象YMCPBに対するSMI制御系を設計するプログラムを次に示します。

%cYMCPB3_smi.m
%-----
 clear all, close all
 dataset=[...
 "ABC_No1_20220128",...
 "ABC_No2_20220131",...
 "ABC_No3_20220128",...
 "ABC_No4_20220131",...
 "ABC_No5_20220131"];
 no=1
 load(dataset(no))
 ID=[8,9,10,11,12,13,14,...
 36,37,38,39,40,41,42,43,...
 66,67,68,69,70,71,72,73,...
 97,98,99,100,101,102,103,...
 126,127,128,129,130,131,132,133];
%-----６次元モデル
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 sys=ss(A0(:,:,ID),B0(:,:,ID),eye(6,6),[]);
%-----３次元モデル
%(1)th,(2)dth,(3)the
 k=[3,6]; j=[2];
 for id=ID
   AA(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,3)];
 end 
 omegae=0.0380*3;
 B=[zeros(2,1);omegae];
 E=eye(3,3);
 C=E(3,:);
 sys3=ss(AA(:,:,ID),B,eye(3,3),[]);  
%-----設計ポイント id=128,ID=33
 Vs=dx_eq_Series(128);    %*3.6
 zs=z_eq_Series(128);     %*100
 ths=th_eq_Series(128);   %/pi*180
 thes=the_eq_Series(128); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,33));
 [A,B]=ssdata(sys3(:,:,33));
%----- 
 AE=[0 -C ;
    zeros(3,1) A];
 BE=[0; B];
 lambda=[-4.7/20 4.7 -2/10];
 nocomp=1;
 specpos=rand(4,3);
 specent=rand(4,3);
 S=swfvpl(AE,BE,lambda,nocomp,specpos,specent);    
 FE=(S*BE)\(S*AE);
 eig(AE-BE*FE);
%---
 [nn,mm]=size(B);
 S1=S(:,1:nn-mm);
 S2=S(:,nn-mm+1:nn);
 Br=eye(nn-mm,mm);
 B2=B(nn-mm+1:nn,:); 
 Lambda=S2*B2;
 Phi=-0.5
 P2=0.5*inv(-Phi);
 Leq=inv(S*BE)*S*AE;
 LPhi=-inv(S*BE)*Phi*S;
 rho=2
 Sr=0
%-----
 L=Leq+LPhi
 Ln=inv(S*BE)*rho
 Lr=inv(Lambda)*(Phi*Sr+S1*Br)
 Ldr=inv(Lambda)*Sr
%-----
 us=thes;
 xs=[0;zs;ths;Vs;0;0];
 x0=[0;0;-1/180*pi;0;0;0];
 yc=0; 
 E=eye(6,6);
 CM=E([3,6],:);
 CM1=100*E([2],:);
 CM2=180/pi*E([3],:);
 ns=2;
%-----
%eof
% L =
%    -0.6402   73.9799   24.2004   11.2935
% Ln =
%    54.0449
% Lr =
%     5.5635
% Ldr =
%      0

図4　３次元モデルに基づくSMI制御系と符号関数を通したSMI制御

ここでは、次の４種類のLQI制御系設計について述べます。

Note e12に示すように、データセットNo.1について、id=128の場合を設計用モデルとして採用したとき、Case2とCase3のどちらを用いても制御性能に大きな差異が見られません。そこでこのコントローラがどれくらいのモデル変動をカバーできるかという意味のロバスト性をシミュレーションにより調べてみます。ここで、モデル変動としては、前進速度を変える場合と、重心位置を変える場合の２通りを考えます。そのために次の３組のデータセットを考えます。

設計用モデルとして、データセットNo.3における船外機取付角-2[deg]のときのid=128の場合を採用し、Case 2とCase 3のLQIコントローラを設計します。そしてデータセットNo.3,4,5における船外機取付角-2[deg]のときのモデル（id=126～133）を安定化できるかを調べます。シミュレーション結果を次に示します。

⇒Case2:いくつかの動作点で不安定
⇒Case3:いくつかの動作点で不安定
図1　データセットNo.3の場合、LQIコントローラのロバスト性

⇒Case2:いくつかの動作点で不安定
⇒Case3:いくつかの動作点で不安定
図2　データセットNo.4の場合、LQIコントローラのロバスト性

⇒Case2:いくつかの動作点で不安定
⇒Case3:いくつかの動作点で不安定
図3　データセットNo.5の場合、LQIコントローラのロバスト性

これらのシミュレーションは本来は非線形シミュレータ上で行うべきですが、ここでは６次元線形モデルを用いて簡易的に行っています。

第１の注意点は、６次元線形モデルの第４番目の状態変数の係数を強制的にとしていることです。これはとなる場合があって、が発散してしまうからです。

第２の注意点は、上のシミュレーションでは各データセットのすべてのモデルを安定化できていませんが、設計用モデルとして、データセットNo.1における船外機取付角-2[deg]のときのid=128の場合を採用すると、各データセットのすべてのモデルは安定化されることです。

これらの事情は、線形モデルが暫定版の非線形シミュレータから得られているので正確ではないためとも考えられます。いずれにしてもできるだけ実機のダイナミックスを反映した非線形シミュレータ上での再検討が必要になります。

%cYMCPB5_robust_lqi.m
%-----
 clear all, close all
 dataset=[...
 "ABC_No1_20220128",...
 "ABC_No2_20220131",...
 "ABC_No3_20220128",...
 "ABC_No4_20220131",...
 "ABC_No5_20220131"];
%-----設計用データセット
 no=3
 load(dataset(no))
 ID=[7,8,9,10,11,12,13,14,...
 36,37,38,39,40,41,42,43,...
 66,67,68,69,70,71,72,73,...
 96,97,98,99,100,101,102,103,...
 126,127,128,129,130,131,132,133];
%-----６次元モデル
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 sys=ss(A0(:,:,ID),B0(:,:,ID),eye(6,6),[]);
%-----５次元モデル
%(1)z,(2)th,(3)dz,(4)dth,(5)the
 k=[2,3,5,6]; i=[2,4,5]; j=[2];
 for id=ID
   A(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,5)];
 end 
 omegae=0.0380*3;
 B=[zeros(4,1);omegae];
 E=eye(5,5);
 CM=E(i,:);
 C=CM(3,:);
 sys5=ss(A(:,:,ID),B,eye(5,5),[]);  
%-----設計ポイント 
 index=35 %[3,11,19,27,35]
 id=ID(index); 
 Vs=dx_eq_Series(id);    %*3.6
 zs=z_eq_Series(id);     %*100
 ths=th_eq_Series(id);   %/pi*180
 thes=the_eq_Series(id); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,index));
 [A,B]=ssdata(sys5(:,:,index));
%-----
 [n,m]=size(B);
 AE=[A B;zeros(m,n+m)];
 BE=[zeros(n,m);eye(m)];
 S=[A B;C 0];
 Mz=0.05;
 Mth=3/180*pi;
 Mthe=10/180*pi; 
 Mue=1;
 Tue=0.05;
 CE=eye(n+m,n+m); 
 QE=diag([1/Mz^2,1/Mth^2,0,0,1/Mthe^2,1/Mue^2]);
 RE=(Tue/Mue)^2;
 [KE,PE]=opt(AE,BE,CE,QE,RE);
 poles=eig(AE-BE*KE);
 FE=KE/S;
%----- 
 F=FE(:,1:n)
 FI=FE(:,n+1:n+m)
 F(1)=0; F(3)=0;
%-----全モデルに対するシミュレーション
for no=[3,4,5] 
for iSIM=index-2:index+5  
 id=ID(iSIM);
 Vs=dx_eq_Series(id);    %*3.6
 zs=z_eq_Series(id);     %*100
 ths=th_eq_Series(id);   %/pi*180
 thes=the_eq_Series(id); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,iSIM));
 Asys(4,4)=0;
 [A,B]=ssdata(sys5(:,:,iSIM));  
 us=thes;
 xs=[0;zs;ths;Vs;0;0];
 x0=[0;0;-1/180*pi;0;0;0];
 yc=0; 
 E=eye(6,6);
 CM=E([2,3,5,6],:);
 CM1=100*E([2],:);
 CM2=180/pi*E([3],:);
 ns=4;
 simout=sim('nYMCPB_lqi.slx');
%[no,iSIM,id,180/pi*us]
 figure(no)
 subplot(221),plot(simout.x.Time,100*simout.x.Data(:,2)),hold on
 subplot(223),plot(simout.y.Time,180/pi*simout.y.Data(:,1)),hold on
 subplot(222),plot(simout.y.Time,180/pi*simout.y.Data(:,3)),hold on
 subplot(224),plot(simout.u.Time,simout.u.Data),hold on
end %iSIM
 figure(no)
 subplot(221),title("z[cm]"),grid
 subplot(223),title("th[deg]"),grid
 subplot(222),title("the[deg]"),grid
 subplot(224),title("ue"),grid
end %no 
%-----
%eof

%cYMCPB3_robust_lqi.m
%-----
 clear all, close all
 dataset=[...
 "ABC_No1_20220128",...
 "ABC_No2_20220131",...
 "ABC_No3_20220128",...
 "ABC_No4_20220131",...
 "ABC_No5_20220131"];
%-----設計用データセット
 no=3
 load(dataset(no))
 ID=[7,8,9,10,11,12,13,14,...
 36,37,38,39,40,41,42,43,...
 66,67,68,69,70,71,72,73,...
 96,97,98,99,100,101,102,103,...
 126,127,128,129,130,131,132,133];
%-----６次元モデル
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 sys=ss(A0(:,:,ID),B0(:,:,ID),eye(6,6),[]);
%-----３次元モデル
%(1)th,(2)dth,(3)the
 k=[3,6]; j=[2];
 for id=ID
   AA(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,3)];
 end 
 omegae=0.0380*3;
 B=[zeros(2,1);omegae];
 E=eye(3,3);
 C=E(3,:);
 sys3=ss(AA(:,:,ID),B,eye(3,3),[]);  
%-----設計ポイント 
 index=35 %[3,11,19,27,35]
 id=ID(index); 
 Vs=dx_eq_Series(id);    %*3.6
 zs=z_eq_Series(id);     %*100
 ths=th_eq_Series(id);   %/pi*180
 thes=the_eq_Series(id); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,index));
 [A,B]=ssdata(sys3(:,:,index));
%-----
 [n,m]=size(B);
 AE=[A B;zeros(m,n+m)];
 BE=[zeros(n,m);eye(m)];
 S=[A B;C 0];
 Mth=3/180*pi;
 Tth=0.5;
 Mthe=5/180*pi; 
 Mue=1;
 Tue=0.05;
 CE=eye(n+m,n+m); 
 QE=diag([1/Mth^2,(Tth/Mth)^2,1/Mthe^2,1/Mue^2]);
 RE=(Tue/Mue)^2;
 [KE,PE]=opt(AE,BE,CE,QE,RE);
 poles=eig(AE-BE*KE);
 FE=KE/S;
%-----  
 F=FE(:,1:n)
 FI=FE(:,n+1:n+m) 
%-----全モデルに対するシミュレーション
for no=[3,4,5] 
for iSIM=index-2:index+5 
 id=ID(iSIM);
 Vs=dx_eq_Series(id);    %*3.6
 zs=z_eq_Series(id);     %*100
 ths=th_eq_Series(id);   %/pi*180
 thes=the_eq_Series(id); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,iSIM));
 Asys(4,4)=0;
 [A,B]=ssdata(sys3(:,:,iSIM));    
 us=thes;
 xs=[0;zs;ths;Vs;0;0];
 x0=[0;0;-1/180*pi;0;0;0];
 yc=0; 
 E=eye(6,6);
 CM=E([3,6],:);
 CM1=100*E([2],:);
 CM2=180/pi*E([3],:);
 ns=2;
 simout=sim('nYMCPB_lqi.slx');
%[no,iSIM,id,180/pi*us]
 figure(no)
 subplot(221),plot(simout.x.Time,100*simout.x.Data(:,2)),hold on
 subplot(223),plot(simout.y.Time,180/pi*simout.y.Data(:,1)),hold on
 subplot(222),plot(simout.y.Time,180/pi*simout.y.Data(:,3)),hold on
 subplot(224),plot(simout.u.Time,simout.u.Data),hold on
end %iSIM
 figure(no)
 subplot(221),title("z[cm]"),grid
 subplot(223),title("th[deg]"),grid
 subplot(222),title("the[deg]"),grid
 subplot(224),title("ue"),grid
end %no 
%-----
%eof

Note e12　LQI制御系の設計
Case 1：５次元モデルに基づくLQI制御

●LQI制御系を設計するために、次の５次元モデルを考えます。

ここで操作入力はではなくとしていますが、制御系の評価の際に

とします。

●制御対象YMCPBに対する制御目的は、あるスラストと取付角の下で、巡航速度で走行するとき、一定の姿勢を保つこととします。制御目的を達成する制御則として次式を考えます。

ここで、右辺第１項は状態フィードバック、第２項は船外機取付角をある収納状態から所定の平衡入力値に設定するための積分動作です。いま

となって、船外機取付角を設定値に維持できたとすると、次式が成り立つ必要があります。

すなわち

ここで、システム行列は正則ですから、任意の平衡入力値に対して、をもたらす状態変数とが定まることが分かります。そこで次の偏差系E3を考えます。

これに対する状態フィードバック

を、二次形式評価関数

を最小化して求めます。(6)は

に注意して

と書けるので、これを積分して次のLQI制御則を得ます。

●制御対象YMCPBに対するLQI制御系を設計するプログラムを次に示します。

%cYMCPB5_lqi.m
%-----
 clear all, close all
 dataset=[...
 "ABC_No1_20220128",...
 "ABC_No2_20220131",...
 "ABC_No3_20220128",...
 "ABC_No4_20220131",...
 "ABC_No5_20220131"];
 no=1
 load(dataset(no))
 ID=[8,9,10,11,12,13,14,...
 36,37,38,39,40,41,42,43,...
 66,67,68,69,70,71,72,73,...
 97,98,99,100,101,102,103,...
 126,127,128,129,130,131,132,133];
%-----６次元モデル
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 sys=ss(A0(:,:,ID),B0(:,:,ID),eye(6,6),[]);
%-----５次元モデル
%(1)z,(2)th,(3)dz,(4)dth,(5)the
 k=[2,3,5,6]; i=[2,4,5]; j=[2];
 for id=ID
   A(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,5)];
 end 
 omegae=0.0380*3;
 B=[zeros(4,1);omegae];
 E=eye(5,5);
 CM=E(i,:);
 C=CM(3,:);
 sys5=ss(A(:,:,ID),B,eye(5,5),[]);  
%-----設計ポイント id=128=ID(33)
 Vs=dx_eq_Series(128);    %*3.6
 zs=z_eq_Series(128);     %*100
 ths=th_eq_Series(128);   %/pi*180
 thes=the_eq_Series(128); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,33));
 [A,B]=ssdata(sys5(:,:,33));
%-----
 [n,m]=size(B);
 AE=[A B;zeros(m,n+m)];
 BE=[zeros(n,m);eye(m)];
 S=[A B;C 0];
 Mz=0.05;
 Mth=3/180*pi;
 Mthe=10/180*pi; 
 Mue=1;
 Tue=0.05;
 CE=eye(n+m,n+m); 
 QE=diag([1/Mz^2,1/Mth^2,0,0,1/Mthe^2,1/Mue^2]);
 RE=(Tue/Mue)^2;
 [KE,PE]=opt(AE,BE,CE,QE,RE);
 poles=eig(AE-BE*KE);
 FE=KE/S;
%----- 
 F=FE(:,1:n)
 FI=FE(:,n+1:n+m)
%F(1)=0; F(3)=0;
%-----  
 us=thes;
 xs=[0;zs;ths;Vs;0;0];
 x0=[0;0;-1/180*pi;0;0;0];
 yc=0; 
 E=eye(6,6);
 CM=E([2,3,5,6],:);
 CM1=100*E([2],:);
 CM2=180/pi*E([3],:);
 ns=4;
%-----
%eof
% F = 0.0728   20.4471    0.4015  -39.3709  186.9394
% FI = 127.8489

図1　LQI制御系シミュレーション

これを実行して次のシミュレーション結果を得ます。

図1　５次元モデルに基づくLQI制御系と符号関数を通したLQI制御

Case 2：Case1でヒーブゲインを零とする場合

図2　ヒーブゲインを零とした場合のLQI制御系と符号関数を通したLQI制御

Case 3：３次元モデルに基づくLQI制御

●制御対象YMCPBに対するLQI制御系を設計するプログラムを次に示します。

%cYMCPB3_lqi.m
%-----
 clear all, close all
 dataset=[...
 "ABC_No1_20220128",...
 "ABC_No2_20220131",...
 "ABC_No3_20220128",...
 "ABC_No4_20220131",...
 "ABC_No5_20220131"];
 no=1
 load(dataset(no))
 ID=[8,9,10,11,12,13,14,...
 36,37,38,39,40,41,42,43,...
 66,67,68,69,70,71,72,73,...
 97,98,99,100,101,102,103,...
 126,127,128,129,130,131,132,133];
%-----６次元モデル
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 sys=ss(A0(:,:,ID),B0(:,:,ID),eye(6,6),[]);
%-----３次元モデル
%(1)th,(2)dth,(3)the
 k=[3,6]; j=[2];
 for id=ID
   AA(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,3)];
 end 
 omegae=0.0380*3;
 B=[zeros(2,1);omegae];
 E=eye(3,3);
 C=E(3,:);
 sys3=ss(AA(:,:,ID),B,eye(3,3),[]);  
%-----設計ポイント id=128=ID(33)
 Vs=dx_eq_Series(128);    %*3.6
 zs=z_eq_Series(128);     %*100
 ths=th_eq_Series(128);   %/pi*180
 thes=the_eq_Series(128); %/pi*180
 [Asys,Bsys,Csys,Dsys]=ssdata(sys(:,:,33));
 [A,B]=ssdata(sys3(:,:,33));
%-----
 [n,m]=size(B);
 AE=[A B;zeros(m,n+m)];
 BE=[zeros(n,m);eye(m)];
 S=[A B;C 0];
 Mth=3/180*pi;
 Tth=0.5
 Mthe=5/180*pi; 
 Mue=1;
 Tue=0.05;
 CE=eye(n+m,n+m); 
 QE=diag([1/Mth^2,(Tth/Mth)^2,1/Mthe^2,1/Mue^2]);
 RE=(Tue/Mue)^2;
 [KE,PE]=opt(AE,BE,CE,QE,RE);
 poles=eig(AE-BE*KE);
 FE=KE/S;
%-----  
 F=FE(:,1:n)
 FI=FE(:,n+1:n+m)
%-----  
 us=thes;
 xs=[0;zs;ths;Vs;0;0];
 x0=[0;0;-1/180*pi;0;0;0];
 yc=0; 
 E=eye(6,6);
 CM=E([3,6],:);
 CM1=100*E([2],:);
 CM2=180/pi*E([3],:);
 ns=2;
%-----
%eof
% F =
%    -2.8411  -45.8237  192.5801
% FI =
%   231.8536

図3　３次元モデルに基づくLQI制御系と符号関数を通したLQI制御

Case 4：Case 1で状態オブザーバを用いる場合

いま、観測変数が状態変数の最初にくるように、状態方程式と観測方程式における状態変数の順番を変えておきます。

このとき、を

のように選んで、オブザーバのパラメータを次のように求めることができます。

ここで、設計パラメータはサイズの行列で

を安定行列とするように選びます。そのためには、仮想システムに対する状態フィードバックによる安定化問題

すなわち

を解けばよいと言えます（、）。

YMCPBの線形モデル

●制御対象YMCPBの運動方程式として次式を考えます。

すなわち

ここで、はそれぞれサージ、ヒーブ、ピッチを、はそれぞれスラスト、船外機取付角を表しています。その他の物理パラメータの説明はここでは省略します。

●いま状態変数ベクトルと操作変数ベクトルをそれぞれ

ととると、次の非線形状態方程式を得ます。

ただし

●制御対象YMCPBに対する制御目的は、あるスラストと取付角の下で、巡航速度で走行するとき、一定の姿勢を保つこととします。これは、次の平衡状態を、次の平衡入力によって保持することとみなすことができます。

この平衡状態と平衡入力を定めるためには、とを所与として、残りのを次の非線形連立方程式を解いて求めます。

したがって、次式が成り立つことに注意します。

●この平衡状態と平衡入力の回りで、非線形状態方程式を線形近似します。

ここで(6)に注意して、線形状態方程式

すなわち

を得ます。

●次の５組のデータセットが与えられています。

各データセットは、133個（）の平衡状態と平衡入力と行列を含みます（行列と行列は固定）。これらは、たとえば、データセットNo1に対しては図1のようにグラフ化されます。

図1 データセットNo1における平衡状態

ここで、赤い〇の平衡状態は、制御系の設計および評価のために選ばれた、次表に示す39点です。の場合の平衡状態回りの線形モデルについてコントローラを設計します。

	run	k
-4 deg	run1	8,9,10,11,12,13,14
0 deg	run1	36,37,38,39,40,41,42,43
4 deg	run1	66,67,68,69,70,71,72,73
8 deg	run1	96,97,98,99,100,101,102,103
2 deg	run1	126,127,128,129,130,131,132,133

またデータセットNo1に対する行列の固有値は図2のようにグラフ化されます。

図2 データセットNo1における行列の固有値

各線形モデルに対して複素固有値が２対、実固有値が２個（うち１個は零）あることがわかります。不安定な複素固有値はピッチの振舞いに、安定な複素固有値はヒーブの振舞いに、非零実固有値は速度変動に関係しています。速度変動が発散する場合があることが気になるところです。

制御系設計用線形モデル

●制御系設計用線形モデルを得るために、まずアクチュエータについては、船外機の取付角を操作するためにレバーを引いた期間だけ一定の角速度で回転するものとします。これを次式で表します。

ここで、は操作入力が正のとき、負のときの値をとるものとします。すなわち

したがって、制御対象の状態方程式(9)とアクチュエータ(10)を合わせて、次式を得ます。

また状態変数のうち、、、が計測できるものとします。このとき、状態方程式(12)に対する観測方程式は次式で表されます。

●制御対象YMCPBに対する制御目的は、あるスラストと取付角の下で、巡航速度で走行するとき、一定の姿勢を保つこととします。この制御目的に照らして、サージ方向の状態変数は省いてよく、また操作変数も船外機取付角に限定することができます。ヒーブ方向の状態変数を残すか省くかによって、次の２種類の線形モデルが考えられます。

５次元モデル

データセットNo1に対する行列の固有値は図3のようにグラフ化されます。ここで、安定な固有値はヒーブの振舞いに、不安定な固有値はピッチの振舞いに関係しています。

図3 データセットNo1を用いた５次元モデルの行列の固有値

３次元モデル

データセットNo1に対する行列の固有値は図4のようにグラフ化されます。これからピッチに対応する固有値は不安定であることが分かります。

図4 データセットNo1を用いた３次元モデルの行列の固有値

●制御目的を達成する制御系を設計する立場から、２つのモデルを比較してみます。まず制御則として次式を考えます。

ここで、右辺第１項は状態フィードバック、第２項は船外機取付角をある収納状態から所定の平衡入力値に設定するための積分動作です。

５次元モデル(14.1)の場合は、ヒーブに関する状態変数を計測できないので、状態オブザーバを用いてこれらを推定して、状態フィードバックを実施することになります。

一方３次元モデル(15.1)の場合は、状態フィードバックは実施できますが、ヒーブの振舞いの影響が外乱として現れますので、これを抑制する必要があります。そのヒーブの振舞いは次式で表されます。

これからヒーブの振舞いには操作入力は直接影響を及ぼさず、ピッチを通して間接的に影響すること分かります。ただ、これは漸近安定（は安定行列）なので、ピッチの安定化を速やかにできれば（が零に収束すれば）、(17)は考慮しなくともよいことになります。しかしながら、３次元モデルに対する状態フィードバックを５次元モデルに適用して、閉ループ系が漸近安定かどうかを調べておくことが考えられます。

制御系設計に３次元モデルを用いる場合は、状態オブザーバが不要となる利点がありますが、ヒーブの振舞いの影響が抑制されていることを確かめるために、５次元モデルに対する閉ループシミュレ－ションを行うことが重要になります。

●以上の結果を、の場合の平衡状態回りの線形モデルについて数値で確認しておきます。

まず元の線形モデル(12)の行列とその固有値は次のようになります。

５次元モデル(14.1)における行列とその固有値は次のようになります。

これからサージに関わる状態変数を除いた場合、ヒーブとピッチに関わる固有値がほぼ引き継がれることが分かります。

一方３次元モデル(15.1)における行列とその固有値は次のようになります。

これからさらにヒーブに関わる状態変数を除いた場合、ピッチに関わる固有値がほぼ引き継がれることが分かります。

●ヒーブとピッチに関わる固有値の可制御性と可観測性を、５次元モデルを用いて調べてみます。

(21)から、すべての固有値について行列の最小特異値は正ですから、形式的には可制御性は成り立つと判定できます。ただ、＊印の値が示すように、ヒーブはピッチに比べて相対的に可制御性の程度が弱く、状態フィードバックによる固有値の変更が難しいと言えます。したがって、制御対象は可制御でなく、むしろ可安定と判定することが考えられます。

一方、(22)から、形式的には可観測性は成り立つと判定できます。ただ、＊印の値が示すように、ヒーブはピッチに比べて相対的に可観測性の程度が弱く、状態オブザーバによる推定が難しいと言えます。

●操作入力（取付角）により、まずピッチが制御され、それがヒーブに影響を与えるという連成の仕組みを調べるために、(15.1)と(17)を数値で確かめてみます。

まず(15.1)は

となり、ヒーブのピッチへの影響を表すに関わる係数の値がそれほど大きくなく、３次元モデルにおいても固有値が継承されることに注意します。

一方(17)は

となります。これは漸近安定なので、ピッチの安定化が速やかに行えれば問題ないとと言えます。

●最後に、取付角からの伝達特性を調べておきます。

図5 取付角からヒーブまで（左図）とピッチまで（右図）の伝達特性

左図は２つの共振特性をもちますが、右図は（２つ目がほとんど消えていて）１つだけの共振特性をもちます。これから、ヒーブのピッチへの連成影響はほとんどないことがわかります。

●以上の数値例の結果を得るためのプログラムを以下に示します。

%ymcpb_ana
%-----
 clear all, close all
 dataset=[...
 "ABC_No1_20220128",...
 "ABC_No2_20220131",...
 "ABC_No3_20220128",...
 "ABC_No4_20220131",...
 "ABC_No5_20220131"];
 no=1
 load(dataset(no))
 ID=[8,9,10,11,12,13,14,...
 36,37,38,39,40,41,42,43,...
 66,67,68,69,70,71,72,73,...
 97,98,99,100,101,102,103,...
 126,127,128,129,130,131,132,133];
 id=1:length(ID);
%-----
 figure(1)
 subplot(321),plot(run_eq_Series,"o-"), hold on
 subplot(321),plot(ID,run_eq_Series(ID),"or"),title("No"+num2str(no)+": "+"run")
 subplot(323),plot(No_eq_Series,"o-"), hold on
 subplot(323),plot(ID,No_eq_Series(ID),"or"),title("No"+num2str(no)+": "+"no")
 subplot(325),plot(dx_eq_Series*3.6,"o-"), hold on
 subplot(325),plot(ID,dx_eq_Series(ID)*3.6,"or"),title("No"+num2str(no)+": "+"dxs")
 subplot(322),plot(z_eq_Series*100,"o-"), hold on
 subplot(322),plot(ID,z_eq_Series(ID)*100,"or"),title("No"+num2str(no)+": "+"zs")
 subplot(324),plot(th_eq_Series/pi*180,"o-"), hold on
 subplot(324),plot(ID,th_eq_Series(ID)/pi*180,"or"),title("No"+num2str(no)+": "+"ths")
 subplot(326),plot(the_eq_Series/pi*180,"o-"), hold on
 subplot(326),plot(ID,the_eq_Series(ID)/pi*180,"or"),title("No"+num2str(no)+": "+"thes")
%-----
 figure(2)
 for i=id, a41(i)=A0(4,1,ID(i)); end
 subplot(331),plot(id,a41,"o-",33,a41(33),"*r"),title("No"+num2str(no)+": "+"a41")
 for i=id, a42(i)=A0(4,2,ID(i)); end
 subplot(332),plot(id,a42,"o-",33,a42(33),"*r"),title("No"+num2str(no)+": "+"a42")
 for i=id, a43(i)=A0(4,3,ID(i)); end
 subplot(333),plot(id,a43,"o-",33,a43(33),"*r"),title("No"+num2str(no)+": "+"a43")
 for i=id, a51(i)=A0(5,1,ID(i)); end
 subplot(334),plot(id,a51,"o-",33,a51(33),"*r"),title("No"+num2str(no)+": "+"a51")
 for i=id, a52(i)=A0(5,2,ID(i)); end
 subplot(335),plot(id,a52,"o-",33,a52(33),"*r"),title("No"+num2str(no)+": "+"a52")
 for i=id, a53(i)=A0(5,3,ID(i)); end
 subplot(336),plot(id,a53,"o-",33,a53(33),"*r"),title("No"+num2str(no)+": "+"a53")
 for i=id, a61(i)=A0(6,1,ID(i)); end
 subplot(337),plot(id,a61,"o-",33,a61(33),"*r"),title("No"+num2str(no)+": "+"a61")
 for i=id, a62(i)=A0(6,2,ID(i)); end
 subplot(338),plot(id,a62,"o-",33,a62(33),"*r"),title("No"+num2str(no)+": "+"a62")
 for i=id, a63(i)=A0(6,3,ID(i)); end
 subplot(339),plot(id,a63,"o-",33,a63(33),"*r"),title("No"+num2str(no)+": "+"a63")
%-----
 figure(3)
 for i=id, a44(i)=A0(4,4,ID(i)); end
 subplot(331),plot(id,a44,"o-",33,a44(33),"*r"),title("No"+num2str(no)+": "+"a44")
 for i=id, a45(i)=A0(4,5,ID(i)); end
 subplot(332),plot(id,a45,"o-",33,a45(33),"*r"),title("No"+num2str(no)+": "+"a45")
 for i=id, a46(i)=A0(4,6,ID(i)); end
 subplot(333),plot(id,a46,"o-",33,a46(33),"*r"),title("No"+num2str(no)+": "+"a46")
 for i=id, a54(i)=A0(5,4,ID(i)); end
 subplot(334),plot(id,a54,"o-",33,a54(33),"*r"),title("No"+num2str(no)+": "+"a54")
 for i=id, a55(i)=A0(5,5,ID(i)); end
 subplot(335),plot(id,a55,"o-",33,a55(33),"*r"),title("No"+num2str(no)+": "+"a55")
 for i=id, a56(i)=A0(5,6,ID(i)); end, axis([0 40 -1 1])
 subplot(336),plot(id,a56,"o-",33,a56(33),"*r"),title("No"+num2str(no)+": "+"a56")
 for i=id, a64(i)=A0(6,4,ID(i)); end, axis([0 40 -1 1])
 subplot(337),plot(id,a64,"o-",33,a64(33),"*r"),title("No"+num2str(no)+": "+"a64")
 for i=id, a65(i)=A0(6,5,ID(i)); end, axis([0 40 -1 1])
 subplot(338),plot(id,a65,"o-",33,a65(33),"*r"),title("No"+num2str(no)+": "+"a65")
 for i=id, a66(i)=A0(6,6,ID(i)); end
 subplot(339),plot(id,a66,"o-",33,a66(33),"*r"),title("No"+num2str(no)+": "+"a66")
%-----
 figure(4)
 for i=id, b41(i)=B0(4,1,ID(i)); end
 subplot(321),plot(id,b41,"o-",33,b41(33),"*r"),title("No"+num2str(no)+": "+"b41")
 for i=id, b42(i)=B0(4,2,ID(i)); end
 subplot(322),plot(id,b42,"o-",33,b42(33),"*r"),title("No"+num2str(no)+": "+"b42")
 for i=id, b51(i)=B0(5,1,ID(i)); end
 subplot(323),plot(id,b51,"o-",33,b51(33),"*r"),title("No"+num2str(no)+": "+"b51")
 for i=id, b52(i)=B0(5,2,ID(i)); end
 subplot(324),plot(id,b52,"o-",33,b52(33),"*r"),title("No"+num2str(no)+": "+"b52")
 for i=id, b61(i)=B0(6,1,ID(i)); end
 subplot(325),plot(id,b61,"o-",33,b61(33),"*r"),title("No"+num2str(no)+": "+"b61")
 for i=id, b62(i)=B0(6,2,ID(i)); end
 subplot(326),plot(id,b62,"o-",33,b62(33),"*r"),title("No"+num2str(no)+": "+"b62")
%-----６次元モデル
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 sys=ss(A0(:,:,ID),B0(:,:,ID),eye(6,6),[]);
 figure(5),pzmap(sys),axis([-10,10,-10,10])
 title("No"+num2str(no)+": "+"k=ID"),hold on
 figure(5),pzmap(sys(:,:,33),"*r")
%-----５次元モデル
%(1)z,(2)th,(3)dz,(4)dth,(5)the
 k=[2,3,5,6]; i=[2,4,5]; j=[2];
 for id=ID
   A(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,5)];
 end 
 B=[zeros(4,1);0.0380*3];
 C=eye(5,5);
 C=C(i,:);
 sys5=ss(A(:,:,ID),B,C,[]);  
 figure(6),pzmap(sys5),axis([-10,10,-10,10])
 title("No"+num2str(no)+": "+"k=ID"),hold on
 figure(6),pzmap(sys5(:,:,33),"*r")
%-----３次元モデル
%(1)th,(2)dth,(3)the
 k=[3,6]; j=[2];
 for id=ID
   AA(:,:,id)=[A0(k,k,id) B0(k,j,id);zeros(1,3)];
 end 
 B=[zeros(2,1);0.0380*3];
 C=eye(3,3);
 sys3=ss(AA(:,:,ID),B,C,[]);  
 figure(7),pzmap(sys3),axis([-10,10,-10,10])
 title("No"+num2str(no)+": "+"k=ID"),hold on
 figure(7),pzmap(sys3(:,:,33),"*r")
%-----設計ポイント id=128,ID=33
 [Asys,Bsys]=ssdata(sys(:,:,33))
 [Asys5,Bsys5,Csys5]=ssdata(sys5(:,:,33))
 [Asys3,Bsys3,Csys3]=ssdata(sys3(:,:,33))
 pl=pole(sys(:,:,33))
 pl5=pole(sys5(:,:,33))
 pl3=pole(sys3(:,:,33))
%-----可制御性 
 w=zeros(5,1);
 for i=1:5
   S=[Bsys5 Asys5-pl5(i)*eye(5,5)];
   w(i)=min(svd(S));
 end
 C0=[pl5 w]
%-----可観測性 
 w=zeros(5,1);
 for i=1:5
   S=[Csys5; Asys5-pl5(i)*eye(5,5)];
   w(i)=min(svd(S));
 end
 OB=[pl5 w]
%-----伝達特性
%(1)x,(2)z,(3)th,(4)dx-V*,(5)dz,(6)dth
 figure(11),hold on
 sigma(sys(:,2,33),logspace(0,2))
 title("No"+num2str(no)+": "+"sys(:,2,33)")
%-----
%(2)z,(5)dz
 figure(12),hold on
 sigma(sys([2,5],2,33),logspace(0,2))
 title("No"+num2str(no)+": "+"sys([2,5],2,33)")
%-----
%(3)th,(6)dth
 figure(13),hold on
 sigma(sys([3,6],2,33),logspace(0,2))
 title("No"+num2str(no)+": "+"sys([3,6],2,33)")
%-----
%eof

参考：漸近安定性

【漸近安定性の定義とその等価な条件】 定義DA:　 条件A1:

漸近安定性とは、定義DAより、平衡状態が乱されたとき復帰できるかどうかにかかわる概念で、条件A1より、行列の固有値がすべて複素左半面にあることを調べて判定されます。

参考：可制御性・可安定性

【可安定性の定義とその等価な条件】 定義DS:　状態フィードバックにより安定化可能 条件S1:　　(はのすべての不安定固有値) 【可制御性の定義とその等価な条件】 定義DC:　任意初期状態を，任意有限時間内に，任意状態に移動可能 条件C2:　 条件C3:　を選んで，の固有値を任意に設定可能 条件C4:　　(はのすべての固有値)

可安定性は、定義DSより、状態フィードバックにより安定化可能であることを意味し、条件S1が判定条件として知られています。一方、可制御性については、条件C4が判定条件として知られています。両者の相違は、すでに左半平面にある安定な固有値について関与（再配置）するかどうかにあります。可安定性の判定では不安定な固有値のみで条件S1を判定しますが、可制御性の判定ではすべて固有値について条件C4を判定します。

それでは、に対して、条件S1と条件C4に出てくる行列

の階数をどのようにして計算するかですが、これはの非零特異値の数で決定します。したがって、条件S1と条件C4は、最小特異値が正かどうかで判定します。

参考：可観測性・可検出性

【可検出性の定義とその等価な条件】 定義DD:　状態オブザーバを構成可能 条件D1:　　(はのすべての不安定固有値) 【可観測性の定義とその等価な条件】 定義DO:　任意有限時間上の入出力データから，初期状態を一意に決定できる 条件O2:　 条件O3:　を選んで，の固有値を任意に設定可能 条件O4:　　(はのすべての固有値)

可検出性は、定義DDより、状態オブザーバを構成可能であることを意味し、条件D1が判定条件として知られています。状態オブザーバの行列はと表され、オブザーバゲインは、仮想システムに対する状態ＦＢによる閉ループ系を安定化して決定します。したがって、可検出性はの可安定性を意味します。一方、可観測性については、条件O4が判定条件として知られています。可検出性との相違は、すでに左半平面にある安定なの固有値についてで再配置するかどうかにあります。可検出性の判定では不安定な固有値のみで条件D1を判定しますが、可観測性の判定ではすべて固有値について条件O4を判定します。

それでは、に対して、条件D1と条件O4に出てくる行列

の階数をどのようにして計算するかですが、これはの非零特異値の数で決定します。したがって、条件D1と条件O4は、最小特異値が正かどうかで判定します。

参考：状態オブザーバ

●通常の状態オブザーバは

で表され

のように漸近的に元の状態を推定するものです。行列はオブザーバゲインと呼ばれ、は安定行列とするように選ばれます。

●最小次元状態オブザーバは

によって表される漸近安定なシステムで、この出力がに漸近するように構成するものです。そのためには、あるサイズの行列に対して次式を満足させることが必要十分となります。

参考：最小実現

一般に，不可制御かつ不可観測な状態空間表現は適当な座標変換により，つぎの正準構造をもつように変換できることが知られています。

ここで，正方行列，，，の次数は一意に定まり，は可制御対，は可観測対です。この正準構造のブロック線図を図1に示します。

図1 正準構造のブロック線図

さて，つぎが成り立ちます。

すなわち，可制御かつ可観測な部分系が入力から出力までの伝達特性を表しています。

これを得るプログラムを、本資料では次のように分割して示しています。

　プログラムIII　************************ 再計画 ************************
・プログラムIII-1　リソース（場所）の定義
・プログラムIII-2　リソース（作業員）の定義
・プログラムIII-3　再計画アクティビティ（仕掛作業、未着手作業）の定義
　　　　　　　　　未着手作業の先行制約の設定
・プログラムIII-4　モードの設定
　　　　　　　　　仕掛作業のモード継承
　　　　　　　　　未着手作業の資源制約（配員）
　　　　　　　　　未着手作業の資源制約（作業場所）
・プログラムIII-5　RCPSP求解

プログラム：Part III

#プログラムIII-1.0
from optseq import *
import math
#=======================================================================
print("")
print(' ************************ リスケ on rdate **********************')
import dill
dill.load_session('usukiR040510b.pkl')  
#=======================================================================
usuki3=Model()
CS=["CS","CC","SS","SC"]

プログラムIII-1

#プログラムIII-1
#=======================================================================
#リソースの定義（場所） 
#=======================================================================
#-----メッシュの使用可否
resPL={}
for id in R:
  resPL[id]=0
  if not id in [4,5,6,13]:
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):        
        resPL[id,i,j] =usuki3.addResource(R[id][0]+"[{0:02d}_{1:02d}]".format(i,j), capacity={(0,"inf"):1}) 
  if id==4:
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in  range(0,11) and j in range(0,12):
          resPL[id,i,j] =usuki3.addResource(R[id][0]+"[{0:02d}_{1:02d}]".format(i,j), capacity={(0,"inf"):0})
        else:     
          resPL[id,i,j] =usuki3.addResource(R[id][0]+"[{0:02d}_{1:02d}]".format(i,j), capacity={(0,"inf"):1})
  if id==5:
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in  range(0,7) and j in range(0,10):
          resPL[id,i,j] =usuki3.addResource(R[id][0]+"[{0:02d}_{1:02d}]".format(i,j), capacity={(0,"inf"):0})
        else:     
          resPL[id,i,j] =usuki3.addResource(R[id][0]+"[{0:02d}_{1:02d}]".format(i,j), capacity={(0,"inf"):1})
  if id==6:       
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in [0,22,23,24,25,26,27,28,29,30,31,32,33,34,35] and j in range(0,13):
          resPL[id,i,j] =usuki3.addResource(R[id][0]+"[{0:02d}_{1:02d}]".format(i,j), capacity={(0,"inf"):0})
        elif i in range(1,17) and j in range(0,1):
          resPL[id,i,j] =usuki3.addResource(R[id][0]+"[{0:02d}_{1:02d}]".format(i,j), capacity={(0,"inf"):0})            
        else:     
          resPL[id,i,j] =usuki3.addResource(R[id][0]+"[{0:02d}_{1:02d}]".format(i,j), capacity={(0,"inf"):1})
  if id==13:       
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in [31,32,33,34,35,36,37,63,64,65,66,67,68,69,70,71] and j in range(0,24):
          resPL[id,i,j] =usuki3.addResource(R[id][0]+"[{0:02d}_{1:02d}]".format(i,j), capacity={(0,"inf"):0})         
        else:     
          resPL[id,i,j] =usuki3.addResource(R[id][0]+"[{0:02d}_{1:02d}]".format(i,j), capacity={(0,"inf"):1})

プログラムIII-2

#プログラムIII-2
#=======================================================================
#リソースの定義（作業員）
#=======================================================================
eps1=1.25
resWH={}    #取付28人
r=math.floor(28*4*eps1);
resWH=usuki3.addResource("R[WH]", capacity={(0,"inf"):r})
#-----
resWW={}    #溶接43人
r=math.floor(43*4*eps1);
resWW=usuki3.addResource("R[WW]", capacity={(0,"inf"):r})
#-----
eps2=1.25
resWF1={}   #鉄艤5人
r=math.floor(5*4*eps2);
resWF1=usuki3.addResource("R[WF1]", capacity={(0,"inf"):r})
#-----
resWF2={}   #甲配7人
r=math.floor(7*4*eps2);
resWF2=usuki3.addResource("R[WF2]", capacity={(0,"inf"):r})
#-----
resWF3={}   #機配5人⇒10人
r=math.floor(10*4*eps2);
resWF3=usuki3.addResource("R[WF3]", capacity={(0,"inf"):r})
#-----
resWP={}    #塗装内4人⇒8人
r=math.floor(8*4*eps2);
resWP=usuki3.addResource("R[WP]", capacity={(0,"inf"):r})
#-----
resWP2={}   #塗装外8人
r=math.floor(8*4*eps2);
resWP2=usuki3.addResource("R[WP2]", capacity={(0,"inf"):r})

プログラムIII-3

#プログラムIII-3.1
#=======================================================================
#再計画 アクティビティ idata3
#=======================================================================
rdate=179
w1=[]; w2=[]; w3=[];
idata12=idata1.union(idata2)
for i in idata12:
  if data[i][3]+data[i][2] < rdate: w1.append(i)
  if data[i][3] <= rdate  and  rdate <= data[i][3]+data[i][2]: w2.append(i)
  if rdate < data[i][3]: w3.append(i)
idata3a=w2;   #仕掛かり中のアクティビティ
idata3b=w3;   #仕掛かっていないアクティビティ 
idata3=w2+w3; #リスケすべきアクティビティ 
print(idata3) 

#プログラムIII-3.2
#=======================================================================
#再計画 作業の定義
#=======================================================================
act3={}
for i in idata3:    
  j=i  
  while data[j][1]!=0:
    j=data[j][4][0]
  act3[i]=usuki3.addActivity(data[i][0], duedate=data[j][3]+data[j][2], \
                                                         backward=True)

#プログラムIII-3.3
#=======================================================================
#再計画 先行制約
#======================================================================= 
for i in idata3b:
  if data[i][1] in ACT_gisou+ACT_kumitate:
    for iw in[0,1,2,3,4]:  
      if data[i][4][iw] in idata3:
        usuki3.addTemporal(act3[i],act3[data[i][4][iw]],\
        tempType=CS[data[i][5][iw][0]], delay= int(data[i][5][iw][1]*2))  
        if data[i][1] in ["A02010","A02610","A03010","A03210","A04310",\
                                     "A09210"] and data[i][5][iw][0]==0:
          usuki3.addTemporal(act3[data[i][4][iw]],act3[i],\
                        tempType=CS[3], delay=-int(data[i][5][iw][1]*2))        
#-----
for i in  idata3b:
  if data[i][1] in PL_gisou+PL_kumitate:
    for iw in [0]:
      if data[i][4][iw] in idata3b:
        ii=data[i][4][iw]      
        usuki3.addTemporal(act3[i],act3[ii],tempType="CS")
        usuki3.addTemporal(act3[ii],act3[i],tempType="SC")       
        for jj in idata3b:   
          for iw in [0]:
            if data[i][4][iw] in idata3b:            
              if data[jj][4][iw]==i:              
                usuki3.addTemporal(act3[jj],act3[i],tempType="CS")
                usuki3.addTemporal(act3[i],act3[jj],tempType="SC")               
#-----開始日の制約
for i in idata3a:   
  usuki3.addTemporal("source",act3[i],tempType="SS",delay= int(data[i][3]))
  usuki3.addTemporal(act3[i],"source",tempType="SS",delay=-int(data[i][3]))     
#-----     
for i in idata3b:
  usuki3.addTemporal("source",act3[i],tempType="SS",delay= rdate)    

プログラムIII-4

#プログラムIII-4.1
#=======================================================================
#仕掛作業 資源制約（配員、場所） 
#=======================================================================
mode3={}
for i in idata3a:  
  if data[i][1] in ACT_gisou:
    mode3[i]=mode1[i] 
    mode3[i].name='mode3['+i+']'
    act3[i].addModes(mode3[i]) 
#-----
  if data[i][1] in [11,1,2,3,4,5,6,12,13,14,15,36,37,38]: #ACT_kumitate.remove(16):  
    mode3[i]=mode2[i] 
    mode3[i].name='mode3['+i+']'    
    act3[i].addModes(mode3[i])
#-----
  if data[i][1] in PL_gisou: 
    mode3[i]=mode1[i] 
    s=mode1[i].name
    mode3[i].name=s[:4]+"3"+s[6:]       
    act3[i].addModes(mode3[i])    
#-----  
  if data[i][1] in PL_kumitate: 
    mode3[i]=mode2[i] 
    s=mode2[i].name
    mode3[i].name=s[:4]+"3"+s[6:]      
    act3[i].addModes(mode3[i])         

#プログラムIII-4.2  
#=======================================================================
#再計画 資源制約（配員）  
#=======================================================================    
for i in idata3b:
  if data[i][1] in ACT_gisou+ACT_kumitate:
#-----
    #1F #2W #3C #4CC #20ブラスト #31鉄艤 #32甲配　#33機配　#34PA #39_甲配 #40_機配   
    if data[i][1] in [1,2,3,4,20,31,32,33,34,39,40]:
      n=math.ceil(data[i][9][0][0]/5)  #１人で何日かかるか
      if data[i][1] in [1,3]:       #取付
        mode3[i]=Mode("mode3["+i+"]",duration=n) 
        mode3[i].addBreak(0,0)
        mode3[i].addResource(resWH,{(0,n):5})
        mode3[i].addParallel(1,n,28)   
      if data[i][1] in [2,4]:       #溶接
        mode3[i]=Mode("mode3["+i+"]",duration=n)    
        mode3[i].addBreak(0,0)
        mode3[i].addResource(resWW,{(0,n):5})
        mode3[i].addParallel(1,n,43)             
      if data[i][1] in [31]:        #鉄艤装
        mode3[i]=Mode("mode3["+i+"]",duration=n)     
        mode3[i].addBreak(0,0)
        mode3[i].addResource(resWF1,{(0,n):5})
        mode3[i].addParallel(1,n,5)          
      if data[i][1] in [32,39]:     #甲板配管
        mode3[i]=Mode("mode3["+i+"]",duration=n)      
        mode3[i].addBreak(0,0)    
        mode3[i].addResource(resWF2,{(0,n):5})
        mode3[i].addParallel(1,n,7)    
      if data[i][1] in [33,40]:     #機関配管
        mode3[i]=Mode("mode3["+i+"]",duration=n)
        mode3[i].addBreak(0,0)    
        mode3[i].addResource(resWF3,{(0,n):5})
        mode3[i].addParallel(1,n,10)     
      if data[i][1]==20:            #塗装内
        mode3[i]=Mode("mode3["+i+"]",duration=n)      
        mode3[i].addBreak(0,0)    
        mode3[i].addResource(resWP, {(0,n):5})   
        mode3[i].addParallel(1,n,8)    
      if data[i][1]==34:            #塗装外
        mode3[i]=Mode("mode3["+i+"]",duration=n)   
        mode3[i].addBreak(0,0)    
        mode3[i].addResource(resWP2,{(0,n):5}) 
        mode3[i].addParallel(1,n,8)    
      act3[i].addModes(mode3[i])
#-----
    #10:入荷 #7:社0 #8:検0 #9:W0 #16:運 #17:運0 #18:運2 #19:合体 #21:運ブ #25:運ブ2	 #35:AT #0:搭載 
    elif data[i][1] in [10,7,8,9,16,17,18,19,21,25,35,0]: 
      mode3[i]=Mode("mode3["+i+"]",duration=data[i][2])   
      act3[i].addModes(mode3[i])
#-----
    #11:開始 #5:社 #6:検 #12:正 #13:反 #14:積み #15:一体 #36:磨き #37:O/P #38:_鉄艤 #16:運 
    elif data[i][1] in [11,5,6,12,13,14,15,36,37,38,16]:      
      mode3[i]=Mode("mode3["+i+"]",duration=data[i][2])  
      act3[i].addModes(mode3[i])       

#プログラムIII-4.3     
#=======================================================================
#再計画 資源制約（作業場所）
#=======================================================================    
for i in idata3b:      
  if data[i][1] in PL_gisou+PL_kumitate:       
#-----
   if data[i][1] in [22,23,24,70,71,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87]: #PL_gisou.remove(72): #72:D5_AH	 
      if data[i][1] in [22,23,24]: data[i][6]=AB
      for no in range(len(data[i][6])):
        no1=data[i][6][no][0];
        L=R[no1][1]; B=R[no1][2];
        L1=data[i][6][no][1]; L2=data[i][6][no][2];
        B1=data[i][6][no][3]; B2=data[i][6][no][4];       
        skipL=math.floor(L/8); skipL=1
        skipB=math.floor(B/2); 
        l=math.floor(data[i][8][0]); b=math.floor(data[i][8][1]);  
        if data[i][1] in [22,23,24]:  l=19; b=16;
        for j in range(L1,L2-l+1,skipL): 
          for k in range(B1,B2-b+1,skipB):
            mode3[i]=Mode("mode3["+i+"]_R[{0:03d}][{1:03d}_{2:03d}][{3:03d}_{4:03d}]".format(no1,j,k,l,b))
            for s in range(0,l):
              for t in range(0,b):
                mode3[i].addBreak(0,0)
                mode3[i].addResource(resPL[no1,j+s,k+t],1,"break")
            act3[i].addModes(mode3[i])            
#-----
    elif data[i][1] in [51,53,54,55,56,57,58,59,60,61,62,63,64]: #PL_kumitate.remove(65): #65:K4_1AS1
      for no in range(len(data[i][6])):
        no1=data[i][6][no][0];
        L=R[no1][1]; B=R[no1][2];
        L1=data[i][6][no][1]; L2=data[i][6][no][2];
        B1=data[i][6][no][3]; B2=data[i][6][no][4];       
        skipL=math.floor(L/8); skipL=1
        skipB=math.floor(B/2); 
        l=math.floor(data[i][8][0]); b=math.floor(data[i][8][1]);  
        if data[i][1] in [22,23,24]:  l=19; b=16;
        for j in range(L1,L2-l+1,skipL): 
          for k in range(B1,B2-b+1,skipB):
            mode3[i]=Mode("mode3["+i+"]_R[{0:03d}][{1:03d}_{2:03d}][{3:03d}_{4:03d}]".format(no1,j,k,l,b))
            for s in range(0,l):
              for t in range(0,b):
                mode3[i].addBreak(0,0)
                mode3[i].addResource(resPL[no1,j+s,k+t],1,"break")
            act3[i].addModes(mode1[i])     
#-----
    elif i in {"A01431","A11701","A11711","A11721"}:  #場所と期間が固定されている場合
      if data[i][6]!=0:
        for no in range(len(data[i][6])):
          no1=data[i][6][no][0];
          L=R[no1][1]; B=R[no1][2];
          L1=data[i][6][no][1]; L2=data[i][6][no][2];
          B1=data[i][6][no][3]; B2=data[i][6][no][4];      
          skipL=math.floor(L/4); skipL=2; 
          skipB=math.floor(B/2); # skipB=1;        
          l=math.floor(data[i][8][0]); b=math.floor(data[i][8][1]);
          for j in range(L1,L2-l+1,skipL):
            for k in range(B1,B2-b+1,skipB):         
              mode3[i]=Mode("mode3["+i+"]_R[{0:03d}][{1:03d}_{2:03d}][{3:03d}_{4:03d}]".format(no1,j,k,l,b),duration=int(data[i][2]))
              for s in range(0,l):
                for t in range(0,b):
                  mode3[i].addResource(resPL[no1,j+s,k+t],1)        
              act3[i].addModes(mode1[i,no1,j,k])                        

プログラムIII-5

#プログラムIII-5
#=======================================================================
#再計画　問題の規模
#======================================================================= 
N=len(usuki3.act)
print("アクティビティ総数（含待機数）:",N)  
M=[]
for a in usuki3.act: M.append(len(a.modes)) 
print("モード数:",M)  
print("平均モード数:",round(sum(M)/N))
P=math.ceil(math.log10(round(sum(M)/N)**N))
print("問題の規模: 10**",P)                       
#=======================================================================
#組立工程計画　問題求解
#=======================================================================
#usuki3.Params.Makespan=True
usuki3.Params.Initial=False
usuki3.Params.TimeLimit=600
usuki3.Params.Neighborhood=20
usuki3.Params.OutputFlag=False
usuki3.optimize()       

プログラムIII-6

#プログラムIII-6
#=======================================================================
#データセットの更新
#=======================================================================
print("")  
for a in usuki3.act:
  for i in idata3b: 
    if a.name==data[i][0]:
      if len(a.modes)==1: nam=a.selected.name
      else: nam=a.selected
      #print("_"+a.name,len(a.modes),nam,a.start,a.completion)
      data[i][3]=a.start
      data[i][2]=a.completion-a.start
      if len(nam)>15:
        r=int(nam[16:19]) 
        x=int(nam[21:24]) 
        y=int(nam[25:28]) 
        l=int(nam[30:33]) 
        b=int(nam[34:37])       
        place=[[ r, x, x+l, y, y+b]]
        data[i][6]=place
        print(data[i][0:8])  
#========================================================================
filename="usu_S1777_resche_R040510c.csv"
print("writing "+filename) 
usuki3.writeExcel(filename)    
import dill
dill.dump_session('usukiR040510c.pkl')  
#------
#eof

組立計画

これを得るプログラムを、本資料では次のように分割して示しています。

　プログラムII　************************ 組立工程計画 ************************
・プログラムII-1　リソース（場所・メッシュ）の定義
・プログラムII-2　リソース（作業員）の定義
・プログラムII-3　組立工程計画アクティビティの定義（一部艤装工程追加）
　　　　　　　　　先行制約の設定
・プログラムII-4　モードの設定
　　　　　　　　　組立工程計画　資源制約（配員）
　　　　　　　　　組立工程計画　資源制約（作業場所）
・プログラムII-5　RCPSP求解
・プログラムII-6　データセットの更新

#プログラムII-1.0
#========================================================================
print("")
print(' ************************ 組立工程計画 *************************')
print("")
#========================================================================
# import dill
# dill.load_session('usukiR040510a.pkl')  
#=====
usuki2=Model()

プログラムII-1

#プログラムII-1.1
#=======================================================================
#リソースの定義（場所） 
#=======================================================================
#-----メッシュの使用可否
for id in R:
  if not id in [4,5,6,13]:
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):        
        resPL[id,i,j] =usuki2.addResource(R[id][0]+"[{0:02d}_{1:02d}]".\
                       format(i,j), capacity=resPL[id,i,j].residual)  
  if id==4:
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in  range(0,11) and j in range(0,12):
          resPL[id,i,j] =usuki2.addResource(R[id][0]+"[{0:02d}_{1:02d}]".\
                         format(i,j), capacity={(0,"inf"):0})
        else:     
          resPL[id,i,j] =usuki2.addResource(R[id][0]+"[{0:02d}_{1:02d}]".\
                         format(i,j), capacity=resPL[id,i,j].residual)  
  if id==5:
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in  range(0,7) and j in range(0,10):
          resPL[id,i,j] =usuki2.addResource(R[id][0]+"[{0:02d}_{1:02d}]".\
                         format(i,j), capacity={(0,"inf"):0})
        else:     
          resPL[id,i,j] =usuki2.addResource(R[id][0]+"[{0:02d}_{1:02d}]".\
                         format(i,j), capacity=resPL[id,i,j].residual)  
  if id==6:       
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in [0,22,23,24,25,26,27,28,29,30,31,32,33,34,35] and \
                                                 j in range(0,13):
          resPL[id,i,j] =usuki2.addResource(R[id][0]+"[{0:02d}_{1:02d}]".\
                         format(i,j), capacity={(0,"inf"):0})
        elif i in range(1,17) and j in range(0,1):
          resPL[id,i,j] =usuki2.addResource(R[id][0]+"[{0:02d}_{1:02d}]".\
                         format(i,j), capacity={(0,"inf"):0})            
        else:     
          resPL[id,i,j] =usuki2.addResource(R[id][0]+"[{0:02d}_{1:02d}]".\
                         format(i,j), capacity=resPL[id,i,j].residual)  
  if id==13:       
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in [31,32,33,34,35,36,37,63,64,65,66,67,68,69,70,71] and \
                                                     j in range(0,24):
          resPL[id,i,j] =usuki2.addResource(R[id][0]+"[{0:02d}_{1:02d}]".\
                         format(i,j), capacity={(0,"inf"):0})         
        else:     
          resPL[id,i,j] =usuki2.addResource(R[id][0]+"[{0:02d}_{1:02d}]".\
                         format(i,j), capacity=resPL[id,i,j].residual)  

プログラムII-2

#プログラムII-2                 
#=======================================================================
#リソースの定義（作業員）
#=======================================================================
resWH =usuki2.addResource("R[WH]",  capacity=resWH.residual)
resWW =usuki2.addResource("R[WW]",  capacity=resWW.residua)
resWF1=usuki2.addResource("R[WF1]", capacity=resWF1.residual)
resWF2=usuki2.addResource("R[WF2]", capacity=resWF2.residual)
resWF3=usuki2.addResource("R[WF3]", capacity=resWF3.residual)
resWP =usuki2.addResource("R[WP]",  capacity=resWP.residual)
resWP2=usuki2.addResource("R[WP2]", capacity=resWP2.residual)

プログラムII-3

#プログラムII-3
#=======================================================================
#組立工程計画アクティビティ idata2
#=======================================================================
w=[]                             
for i in data:
 if data[i][1] in ACT_kumitate+PL_kumitate and data[i][11] in BLK_all and data[i][12] in [0]:
   w.append(i)
idata2=w&data.keys()    
if 3 in BLK_all:
  idata2.add("A09320") #組立中ブラスト（運ブ）
  idata2.add("A09321") #組立中ブラスト（ブラスト棟）
  idata2.add("A09322") #組立中ブラスト（ブラスト作業）
  idata2.add("A09330") #組立中ブラスト（運ブ２）
if 7 in BLK_all:
  idata2=idata2-{
  'A11700','A11701','A11702','A11703','A11704','A11705','A11706',\
  'A11710','A11711','A11712','A11713','A11714','A11715','A11716',\
  'A11720','A11721','A11722','A11723','A11724','A11725','A11726','A11727','A11728','A11729'} #総組は場所だけ残して作業は除く
print(idata2)

プログラムII-4.1

#プログラムII-4.1a
#=======================================================================
#組立工程計画　作業の定義
#=======================================================================
act2={}
for i in idata2:    
  j=i  
  while data[j][1]!=16:
    j=data[j][4][0]
  act2[i]=usuki2.addActivity(data[i][0], duedate=data[j][3]+data[j][2], backward=True)      

#プログラムII-4.1b  
#=======================================================================
#組立工程計画　先行制約
#=======================================================================
for i in idata2:
  if data[i][1] in ACT_kumitate:
    for iw in[0,1,2,3,4]:
      if data[i][4][iw] in idata2:
        usuki2.addTemporal(act2[i],act2[data[i][4][iw]],tempType=CS[data[i][5][iw][0]], delay= int(data[i][5][iw][1]*2))  
        if data[i][1] in ["A02010","A02610","A03010","A03210","A04310","A09210",] and data[i][5][iw][0]==0: #外板
          usuki2.addTemporal(act2[data[i][4][iw]],act2[i],tempType=CS[3], delay=-int(data[i][5][iw][1]*2))   
#-----
for i in idata2:
  if data[i][1] in PL_kumitate:
    for iw in[0]:
      if data[i][4][iw] in idata2:
        ii=data[i][4][iw]      
        usuki2.addTemporal(act2[i],act2[ii],tempType="CS")
        usuki2.addTemporal(act2[ii],act2[i],tempType="SC")       
        for jj in idata2:   
          for iw in[0]:
            if data[i][4][iw] in idata2:            
              if data[jj][4][iw]==i:              
                usuki2.addTemporal(act2[jj],act2[i],tempType="CS")
                usuki2.addTemporal(act2[i],act2[jj],tempType="SC")   
#-----開始日の制約
for i in idata2:
  if data[i][1] in [11]:
    usuki2.addTemporal("source",act2[i],tempType="SS",delay= 80)    
#-----運搬日の制約
for i in idata2:
  if data[i][1] in [16]:
    usuki2.addTemporal("source",act2[i],tempType="SS",delay= int(data[i][3]-10))      

プログラムII-4.2

#プログラムII-4.2a  
#=======================================================================
#組立工程計画 　資源制約（配員）  
#=======================================================================
mode2={}
for i in idata2:
  if data[i][1] in ACT_kumitate+ACT_gisou:
#-----
    #1F #2W #3C #4CC #20ブラスト #31鉄艤 #32甲配　#33機配　#34PA #39_甲配 #40_機配   
    if data[i][1] in [1,2,3,4,20,31,32,33,34,39,40]:
      n=math.ceil(data[i][9][0][0]/5)  #１人で何日かかるか
      if data[i][1] in [1,3]:       #取付
        mode2[i]=Mode("mode2["+i+"]",duration=n) 
        mode2[i].addBreak(0,0)
        mode2[i].addResource(resWH,{(0,n):5})
        mode2[i].addParallel(1,n,28)   
      if data[i][1] in [2,4]:       #溶接
        mode2[i]=Mode("mode2["+i+"]",duration=n)    
        mode2[i].addBreak(0,0)
        mode2[i].addResource(resWW,{(0,n):5})
        mode2[i].addParallel(1,n,43)             
      if data[i][1] in [31]:        #鉄組立
        mode2[i]=Mode("mode0["+i+"]",duration=n)     
        mode2[i].addBreak(0,0)
        mode2[i].addResource(resWF1,{(0,n):5})
        mode2[i].addParallel(1,n,5)          
      if data[i][1] in [32,39]:     #甲板配管
        mode2[i]=Mode("mode0["+i+"]",duration=n)      
        mode2[i].addBreak(0,0)    
        mode2[i].addResource(resWF2,{(0,n):5})
        mode2[i].addParallel(1,n,7)    
      if data[i][1] in [33,40]:     #機関配管
        mode2[i]=Mode("mode0["+i+"]",duration=n)
        mode2[i].addBreak(0,0)    
        mode2[i].addResource(resWF3,{(0,n):5})
        mode2[i].addParallel(1,n,10)     
      if data[i][1]==20:            #塗装内
        mode2[i]=Mode("mode0["+i+"]",duration=n)      
        mode2[i].addBreak(0,0)    
        mode2[i].addResource(resWP, {(0,n):5})   
        mode2[i].addParallel(1,n,8)    
      if data[i][1]==34:            #塗装外
        mode2[i]=Mode("mode0["+i+"]",duration=n)   
        mode2[i].addBreak(0,0)    
        mode2[i].addResource(resWP2,{(0,n):5}) 
        mode2[i].addParallel(1,n,8)    
      act2[i].addModes(mode2[i])
#-----
    #10:入荷 #7:社0 #8:検0 #9:W0 #16:運 #17:運0 #18:運2 #19:合体 #21:運ブ #25:運ブ2 #35:AT #0:搭載 
    elif data[i][1] in [10,7,8,9,16,17,18,19,21,25,35,0]: 
      mode2[i]=Mode("mode2["+i+"]",duration=data[i][2]) #組立工程運搬作業時数未定のため期間だけ確保 
      act2[i].addModes(mode2[i])
#-----
    #11:開始 #5:社 #6:検 #12:正 #13:反 #14:積み #15:一体 #36:磨き #37:O/P #38:_鉄艤 #16:運 
    elif data[i][1] in [11,5,6,12,13,14,15,36,37,38,16]:      
      mode2[i]=Mode("mode2["+i+"]",duration=data[i][2]) #組立工程運搬作業時数未定のため期間だけ確保   
      act2[i].addModes(mode2[i])                

#プログラムII-4.2b                                
#=======================================================================
#組立工程計画 資源制約（作業場所）      
#=======================================================================
for i in idata2:      
  if data[i][1] in PL_kumitate+PL_gisou:   
#-----    
    if data[i][1] in [22,23,24,70,71,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87]: #PL_gisou.remove(72): #72:D5_AH	 
      if data[i][1] in [22,23,24]: data[i][6]=AB
      for no in range(len(data[i][6])):
        no1=data[i][6][no][0];
        L=R[no1][1]; B=R[no1][2];
        L1=data[i][6][no][1]; L2=data[i][6][no][2];
        B1=data[i][6][no][3]; B2=data[i][6][no][4];       
        skipL=math.floor(L/8); skipL=1
        skipB=math.floor(B/2); 
        l=math.floor(data[i][8][0]); b=math.floor(data[i][8][1]);  
        if data[i][1] in [22,23,24]:  l=19; b=16;
        for j in range(L1,L2-l+1,skipL): 
          for k in range(B1,B2-b+1,skipB):
            mode2[i]=Mode("mode2["+i+"]_R[{0:03d}][{1:03d}_{2:03d}][{3:03d}_{4:03d}]".format(no1,j,k,l,b))
            for s in range(0,l):
              for t in range(0,b):
                mode2[i].addBreak(0,0)
                mode2[i].addResource(resPL[no1,j+s,k+t],1,"break")
            act2[i].addModes(mode2[i])            
#-----
    elif data[i][1] in [51,53,54,55,56,57,58,59,60,61,62,63,64]: #PL_kumitate.remove(65): #65:K4_1AS1
      for no in range(len(data[i][6])):
        no1=data[i][6][no][0];
        L=R[no1][1]; B=R[no1][2];
        L1=data[i][6][no][1]; L2=data[i][6][no][2];
        B1=data[i][6][no][3]; B2=data[i][6][no][4];       
        skipL=math.floor(L/8); skipL=1
        skipB=math.floor(B/2); 
        l=math.floor(data[i][8][0]); b=math.floor(data[i][8][1]);  
        if data[i][1] in [22,23,24]:  l=19; b=16;
        for j in range(L1,L2-l+1,skipL): 
          for k in range(B1,B2-b+1,skipB):
            mode2[i,no1,j,k]=Mode("mode2["+i+"]_R[{0:03d}][{1:03d}_{2:03d}][{3:03d}_{4:03d}]".format(no1,j,k,l,b))
            for s in range(0,l):
              for t in range(0,b):
                mode2[i].addBreak(0,0)
                mode2[i].addResource(resPL[no1,j+s,k+t],1,"break")
            act2[i].addModes(mode2[i])     
#-----
    elif i in {"A01431","A11701","A11711","A11721"}:  #場所と期間が固定されている場合
      if data[i][6]!=0:
        for no in range(len(data[i][6])):
          no1=data[i][6][no][0];
          L=R[no1][1]; B=R[no1][2];
          L1=data[i][6][no][1]; L2=data[i][6][no][2];
          B1=data[i][6][no][3]; B2=data[i][6][no][4];      
          skipL=math.floor(L/4); skipL=2; 
          skipB=math.floor(B/2); # skipB=1;        
          l=math.floor(data[i][8][0]); b=math.floor(data[i][8][1]);
          for j in range(L1,L2-l+1,skipL):
            for k in range(B1,B2-b+1,skipB):         
              mode2[i]=Mode("mode2["+i+"]_R[{0:03d}][{1:03d}_{2:03d}][{3:03d}_{4:03d}]".format(no1,j,k,l,b),duration=int(data[i][2]))
              for s in range(0,l):
                for t in range(0,b):
                  mode2[i].addResource(resPL[no1,j+s,k+t],1)        
              act2[i].addModes(mode2[i])                

プログラムII-5

#プログラムII-5
#=======================================================================
#組立工程計画　問題の規模
#=======================================================================
N=len(usuki2.act)
print("アクティビティ総数（含待機数）:",N)  
M=[]
for a in usuki2.act: M.append(len(a.modes)) 
print("モード数:",M)  
print("平均モード数:",round(sum(M)/N))
P=math.ceil(math.log10(round(sum(M)/N)**N))
print("問題の規模: 10**",P)                       
#=======================================================================
#組立工程計画　問題求解
#=======================================================================
#usuki2.Params.Makespan=True
usuki2.Params.Initial=False
usuki2.Params.TimeLimit=3600
usuki2.Params.Neighborhood=20
usuki2.Params.OutputFlag=False
usuki2.optimize()      

●プログラムII-6

#プログラムII-6             
#=======================================================================
#データセットの更新
#=======================================================================
print("")  
for a in usuki2.act:
  for i in idata2: 
    if a.name==data[i][0]:
      if len(a.modes)==1: nam=a.selected.name
      else: nam=a.selected
     #print("_"+a.name,len(a.modes),nam,a.start,a.completion)
      data[i][3]=a.start
      data[i][2]=a.completion-a.start
      if len(nam)>15:
        r=int(nam[16:19]) 
        x=int(nam[21:24]) 
        y=int(nam[25:28]) 
        l=int(nam[30:33]) 
        b=int(nam[34:37])       
        place=[[ r, x, x+l, y, y+b]]
        data[i][6]=place
        print(data[i][0:8])  
#========================================================================
filename="usu_S1777_kumitate_R040510b.csv"
print("writing "+filename) 
usuki2.writeExcel(filename)  
import dill
dill.dump_session('usukiR040510b.pkl')  
#------
#eof

艤装計画

プログラム：Part I-4

●プログラムI-4.1

#プログラムI-4.1
#=======================================================================
#艤装工程計画 　資源制約（配員）  
#=======================================================================
mode1={}
for i in idata1:  
  if data[i][1] in ACT_gisou:
#-----
    #1F #2W #3C #4CC #20ブラスト #31鉄艤 #32甲配　#33機配　#34PA #39_甲配 #40_機配   
    if data[i][1] in [1,2,3,4,20,31,32,33,34,39,40]:
      n=math.ceil(data[i][9][0][0]/5)  #１人で何日かかるか
      if data[i][1] in [1,3]:       #取付
        mode1[i]=Mode("mode0["+i+"]",duration=n) 
        mode1[i].addBreak(0,0)
        mode1[i].addResource(resWH,{(0,n):5})
        mode1[i].addParallel(1,n,28)   
      if data[i][1] in [2,4]:       #溶接
        mode1[i]=Mode("mode0["+i+"]",duration=n)    
        mode1[i].addBreak(0,0)
        mode1[i].addResource(resWW,{(0,n):5})
        mode1[i].addParallel(1,n,43)             
      if data[i][1] in [31]:        #鉄艤装
        mode1[i]=Mode("mode2["+i+"]",duration=n)     
        mode1[i].addBreak(0,0)
        mode1[i].addResource(resWF1,{(0,n):5})
        mode1[i].addParallel(1,n,5)          
      if data[i][1] in [32,39]:     #甲板配管
        mode1[i]=Mode("mode2["+i+"]",duration=n)      
        mode1[i].addBreak(0,0)    
        mode1[i].addResource(resWF2,{(0,n):5})
        mode1[i].addParallel(1,n,7)    
      if data[i][1] in [33,40]:     #機関配管
        mode1[i]=Mode("mode2["+i+"]",duration=n)
        mode1[i].addBreak(0,0)    
        mode1[i].addResource(resWF3,{(0,n):5})
        mode1[i].addParallel(1,n,10)     
      if data[i][1]==20:            #塗装内
        mode1[i]=Mode("mode2["+i+"]",duration=n)      
        mode1[i].addBreak(0,0)    
        mode1[i].addResource(resWP, {(0,n):5})   
        mode1[i].addParallel(1,n,8)    
      if data[i][1]==34:            #塗装外
        mode1[i]=Mode("mode2["+i+"]",duration=n)   
        mode1[i].addBreak(0,0)    
        mode1[i].addResource(resWP2,{(0,n):5}) 
        mode1[i].addParallel(1,n,8)    
      act1[i].addModes(mode1[i])
#-----
    #10:入荷 #7:社0 #8:検0 #9:W0 #16:運 #17:運0 #18:運2 #19:合体 #21:運ブ #25:運ブ2 #35:AT #0:搭載 
    elif data[i][1] in [10,7,8,9,16,17,18,19,21,25,35,0]: 
      mode1[i]=Mode("mode2["+i+"]",duration=data[i][2]) #艤装工程運搬作業時数未定のため期間だけ確保 
      act1[i].addModes(mode1[i])
#-----
    #11:開始 #5:社 #6:検 #12:正 #13:反 #14:積み #15:一体 #36:磨き #37:O/P #38:_鉄艤 #16:運 
    elif data[i][1] in [11,5,6,12,13,14,15,36,37,38,16]:      
      mode1[i]=Mode("mode0["+i+"]",duration=data[i][2]) #組立工程運搬作業時数未定のため期間だけ確保   
      act1[i].addModes(mode1[i])       

●プログラムI-4.2

#プログラムI-4.2
#=======================================================================
#艤装工程計画 資源制約（作業場所）      
#=======================================================================
for i in idata1:      
  if data[i][1] in PL_gisou:      
#----- 
    if data[i][1] in [22,23,24,70,71,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87]: #PL_gisou.remove(72): #72:D5_AH	 
      if data[i][1] in [22,23,24]: data[i][6]=AB
      for no in range(len(data[i][6])):
        no1=data[i][6][no][0];
        L=R[no1][1]; B=R[no1][2];
        L1=data[i][6][no][1]; L2=data[i][6][no][2];
        B1=data[i][6][no][3]; B2=data[i][6][no][4];       
        skipL=math.floor(L/8); skipL=1
        skipB=math.floor(B/2); 
        l=math.floor(data[i][8][0]); b=math.floor(data[i][8][1]);  
        if data[i][1] in [22,23,24]:  l=19; b=16;
        for j in range(L1,L2-l+1,skipL): 
          for k in range(B1,B2-b+1,skipB):
            mode1[i]=Mode("mode2["+i+"]_R[{0:03d}][{1:03d}_{2:03d}][{3:03d}_{4:03d}]".format(no1,j,k,l,b))
            for s in range(0,l):
              for t in range(0,b):
                mode1[i].addBreak(0,0)
                mode1[i].addResource(resPL[no1,j+s,k+t],1,"break")
            act1[i].addModes(mode1[i])            
#-----
    elif data[i][1] in [51,53,54,55,56,57,58,59,60,61,62,63,64]: #PL_kumitate.remove(65): #65:K4_1AS1
      for no in range(len(data[i][6])):
        no1=data[i][6][no][0];
        L=R[no1][1]; B=R[no1][2];
        L1=data[i][6][no][1]; L2=data[i][6][no][2];
        B1=data[i][6][no][3]; B2=data[i][6][no][4];       
        skipL=math.floor(L/8); skipL=1
        skipB=math.floor(B/2); 
        l=math.floor(data[i][8][0]); b=math.floor(data[i][8][1]);  
        if data[i][1] in [22,23,24]:  l=19; b=16;
        for j in range(L1,L2-l+1,skipL): 
          for k in range(B1,B2-b+1,skipB):
            mode1[i]=Mode("mode2["+i+"]_R[{0:03d}][{1:03d}_{2:03d}][{3:03d}_{4:03d}]".format(no1,j,k,l,b))
            for s in range(0,l):
              for t in range(0,b):
                mode1[i].addBreak(0,0)
                mode1[i].addResource(resPL[no1,j+s,k+t],1,"break")
            act1[i].addModes(mode1[i])     
#-----
    elif i in {"A01431","A11701","A11711","A11721"}:  #場所と期間が固定されている場合(1AS1,D5_AH)
      if data[i][6]!=0:
        for no in range(len(data[i][6])):
          no1=data[i][6][no][0];
          L=R[no1][1]; B=R[no1][2];
          L1=data[i][6][no][1]; L2=data[i][6][no][2];
          B1=data[i][6][no][3]; B2=data[i][6][no][4];      
          skipL=math.floor(L/4); skipL=2; 
          skipB=math.floor(B/2); # skipB=1;        
          l=math.floor(data[i][8][0]); b=math.floor(data[i][8][1]);
          for j in range(L1,L2-l+1,skipL):
            for k in range(B1,B2-b+1,skipB):         
              mode1[i]=Mode("mode2["+i+"]_R[{0:03d}][{1:03d}_{2:03d}][{3:03d}_{4:03d}]".format(no1,j,k,l,b),duration=int(data[i][2]))
              for s in range(0,l):
                for t in range(0,b):
                  mode1[i].addResource(resPL[no1,j+s,k+t],1)        
              act1[i].addModes(mode1[i])       

●プログラムI-5

#プログラムI-5  
#=======================================================================
#艤装工程計画 問題の規模
#=======================================================================
N=len(usuki1.act)
print("アクティビティ総数（含待機数）:",N)  
M=[]
for a in usuki1.act: M.append(len(a.modes))
print("モード数:",M)  
print("平均モード数:",round(sum(M)/N))
P=math.ceil(math.log10(round(sum(M)/N)**N))
print("問題の規模: 10**",P)                       
#=======================================================================
#艤装工程計画 問題求解
#=======================================================================
usuki1.Params.Makespan=True
usuki1.Params.Initial=False
usuki1.Params.TimeLimit=3600
usuki1.Params.Neighborhood=20
# usuki1.Params.RandomSeed=10
usuki1.Params.OutputFlag=False
usuki1.optimize()

●プログラムI-6

#プログラムI-6             
#=======================================================================
#データセットの更新
#=======================================================================
print("")  
for a in usuki1.act:
  for i in idata1: 
    if a.name==data[i][0]:
      if len(a.modes)==1: nam=a.selected.name
      else: nam=a.selected
     #print("_"+a.name,len(a.modes),nam,a.start,a.completion)
      data[i][3]=a.start
      data[i][2]=a.completion-a.start
      if len(nam)>15:
        r=int(nam[16:19]) 
        x=int(nam[21:24]) 
        y=int(nam[25:28]) 
        l=int(nam[30:33]) 
        b=int(nam[34:37])       
        place=[[ r, x, x+l, y, y+b]]
        data[i][6]=place
        print(data[i][0:8])  
#========================================================================
filename="usu_S1777_gisou_R040510a.csv"
print("writing "+filename)  
usuki1.writeExcel(filename)   
# import dill
# dill.dump_session('usukiR040510a.pkl')  
# sys.exit()     
#-----
#eof

現行計画

現行計画（総合日程表）における各ブロックのブロック搭載日、ブラスト開始日、組立開始日また外注の場合の入荷日をグラフにしたものを次に示します。

「工程’s」にはリソースの山積機能がありますので、これを利用して作業時数と場所の専有面積を山積した結果を次に示します。

時数山積み

面積山積み

実際の配員計画と定盤計画では様々なやりくり（残業・休日出勤、段積み・先行搭載など）がなされていて、実際の造船工程が回っていると言えます。

RCPSP計画
当該造船所の工程計画問題をRCPSPとして定式化した場合、大規模な組合せ最適化問題となり、ソルバーOptSeqを用いても、問題全体を一括して解くことは困難です。そのため、造船工程をブロック組立工程と艤装・待機工程に分割して、適切につなぎ合わせる工夫が必要となります。一般にブロック搭載工程は固定してよいので、下流の艤装・待機工程はプル型すなわち搭載日程に後詰めします。これに、上流の組立工程をジャストインすることになりますが、その際前詰め（方針I）するか後詰め（方針II）するかの２通りが考えられます。各ブロック組立の開始順番などは経験的に確立されているとしてよいので、総合日程表におけるブロック開始日を固定してこれに前詰めするのが方針Iです。ただ、これはジャストインの達成は難しいので、接続を行う特別の場所の確保が必要となります。これらの方針を次図に示します。

一方、組立工程を艤装・待機工程にジャストインする方針を次図に示します。

どちらともリソース（作業員と作業場所）が分離されていることが前提となりますが、実際には一部例外的な処理が必要です。

また共通した重要な課題として、日程計画・配員計画・定盤計画（場所の割付）をどう連動させるかがあります。本研究では、OptSeqのもつ小作業並列化機能を用いて作業時数と配員制約から作業期間を決めます。同時にその未知の期間に場所の割付ができるかを、OptSeqのもつ仮想アクティビティの手法を用いてチェックします。すなわち、OptSeqと用いて日程計画・配員計画・定盤計画の３つを連動させることを試みます。

●種々検討の結果、上の現行計画に対応して、方針IIに基づいて次図を得ることができました。

艤装工程計画を得るプログラムを、本資料では次のように分割して示しています。

　プログラム I　************************ 艤装工程計画 ************************
・プログラム I-1　リソース（場所・メッシュ、作業員）の定義
・プログラム I-2　データセットの定義
・プログラム I-3　艤装工程計画アクティビティの定義（一部組立工程追加）
　　　　　　　　　先行制約の設定
・プログラム I-4　モードの設定
　　　　　　　　　艤装工程計画　資源制約（配員）
　　　　　　　　　艤装工程計画　資源制約（作業場所）　
・プログラム I-5　RCPSP求解
・プログラム I-6　データセットの更新

プログラム：Part I-3

#プログラムI-3.0
#=======================================================================
#識別番号 
#=======================================================================
#作業項目・作業場所の識別番号 data[i][1]
#----- 
#11:開始 #1:F #2:W #3:C #4:CC #5:社 #6:検 #12:正 #13:反 #14:積み #15:一体 #36:磨き #37:O/P #38:_鉄艤 #16:運 
ACT_kumitate=[11,1,2,3,4,5,6,12,13,14,15,36,37,38,16]
#----- 
#10:入荷 #7:社0 #8:検0 #9:W0 #16:運 #17:運0 #18:運2 #19:合体 #20:ブラスト #21:運ブ #25:運ブ2	
#31:鉄艤 #32:甲配 #33:機配 #39:_甲配 #40:_機配 #34:PA #35:AT #0:搭載 
ACT_gisou=[10,7,8,9,16,17,18,19,20,21,25,31,32,33,39,40,34,35,0] 
#----- 
#51:K1K2 #53:K3	 #54:K4 #55:K5U	 #55:K5R #56:K6	 
#57:K1_1FS1 #58:K1_1S6P #59:K1_1S6S #60:K2_8S6P #61:K2_8S6S #62:K4S4 #63:県岸 #64:K1K2_K5U #65:K4_1AS1
PL_kumitate=[51,53,54,55,56,57,58,59,60,61,62,63,64,65]		
#----- 
#22:AB #23:BB #24:CB #70:RD #71:RD2 #72:D5_AH		
#73:D7_3AS1 #74:D7_3AS12 #75:D7_1AS1P #76:D7_1AS1P2 #77:D7_1AS1S #78:D7_1AS1S2 #79:D7_1AS5	
#80:D7_1AS5P2 #81:D7_1AS5S #82:D7_1AS5S2 #83:D7_2S6P #84:D7_2S6S #85:D7_AH1P #86:D7_AH1S #87:D7_AH
PL_gisou=[22,23,24,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87]
#----- 
#99:その他
ACT_misc=[99]
#----- 
#ブロックカテゴリの識別番号 data[i][11] 
#1:S1 #2:AS #3:FS #4:S6 #5:L1 #6:D1 #7:AH 
BLK_all=[1,2,3,4,5,6,7]
#----- 
#外注の識別番号 data[i][12] 
#1:外注 #2:内製
ORD_all=[1,0]

●プログラムI-3.1

#プログラムI-3.1a
#=======================================================================
#艤装工程計画アクティビティのidata1
#=======================================================================
w=[]                             
for i in data:
  if (data[i][1] in ACT_gisou+PL_gisou) and (data[i][11] in BLK_all) and (data[i][12] in ORD_all):
    w.append(i)
idata1=w&data.keys()
#-----
if 7 in BLK_all:
  idata1=idata1-{\
  'A11700','A11701','A11702','A11703','A11704','A11705','A11706',\
  'A11710','A11711','A11712','A11713','A11714','A11715','A11716',\
  'A11720','A11721','A11722','A11723','A11724','A11725','A11726','A11727','A11728','A11729'} #総組は場所さけ残して作業は除く
print(idata1)
#====
act1={}
for i in idata1:   
    j=i  
    while not (data[j][1] in [0]):
      j=data[j][4][0]
    act1[i]=usuki1.addActivity(data[i][0], duedate=data[j][3]+data[j][2], backward=True)           

●プログラムI-3.1b
#=======================================================================
#艤装工程計画 先行制約
#=======================================================================
for i in  idata1:
  if data[i][1] in ACT_gisou:
    for iw in[0,1,2,3,4]:
      if data[i][4][iw] in idata1:
        usuki1.addTemporal(act1[i],act1[data[i][4][iw]],tempType=CS[data[i][5][iw][0]], delay= int(data[i][5][iw][1]*2))  
        if data[i][1] in ["A02010","A02610","A03010","A03210","A04310","A09210",] and data[i][5][iw][0]==0: #外板
          usuki1.addTemporal(act1[data[i][4][iw]],act1[i],tempType=CS[3], delay=-int(data[i][5][iw][1]*2))   
#-----
for i in  idata1:
  if data[i][1] in PL_gisou:
    for iw in [0]:
      if data[i][4][iw] in idata1:
        ii=data[i][4][iw]      
        usuki1.addTemporal(act1[i],act1[ii],tempType="CS")
        usuki1.addTemporal(act1[ii],act1[i],tempType="SC")       
        for jj in idata1:   
          for iw in [0]:
            if data[i][4][iw] in idata1:            
              if data[jj][4][iw]==i:              
                usuki1.addTemporal(act1[jj],act1[i],tempType="CS")
                usuki1.addTemporal(act1[i],act1[jj],tempType="SC")              
#-----搭載日の固定
for i in idata1:
  if data[i][1] in [0,11]:
    usuki1.addTemporal("source",act1[i],tempType="SS",delay= int(data[i][3]))
    if data[i][1]==0:
      usuki1.addTemporal(act1[i],"source",tempType="SS",delay=-int(data[i][3]))     

S1777ブロック一覧

番船S1777では次の119個のブロックを扱います。

総合日程表の「工程’s」による表現
●RCPSP計画に必要なデータセットの作成のために、まず総合日程表を「工程’s」で表します。そこでは資源マスターとして次表を定義しています。

データセット
●RCPSP計画に必要なデータセットは、総合日程表を「工程’s」で表したあと、「エクスポート」機能を用いてTSVファイルを得て、これを編集して生成します。

●データセットは1847 個のレコードから構成され、１レコードの各フィールドを次に示します。

右列は、次の一例を対応させたものです。

“A00602”:[“S1777_4S1S_A00602_14Lx12B_F”,1,8,83 ,[“A00603″,”A00702″,”0″,”0″,”0”],[[0,0],[ 2,0],[ 0,0],[ 0,0],[ 0,0]], 0, 0,[ 14 , 12 , 73 ],[[152 ],[0 ],[0 ]],[[0 ],[0 ],[0 ]],1, 0, ],\

上記フィールドで作業番号jobは次表で定義されます。

●作業時数については「工数管理システム」においてログが取られています。これに基づいて標準作業時間が決められています。

プログラム：Part I-2

●プログラムI-2.1

#プログラムI-2.1a
#=====データセット
#id":["0:番船_ブロック_id_サイズ_作業",1:job,2:期間,3:開始,
#["4:後続＃1"," 4:後続＃2"," 4:後続＃3"," 4:後続＃4"],
#[[5:型#1,5:型#2],[ 5:型#1,5:型#2],[ 5:型#1,5:型#2],[ 5:型#1,5:型#2]], 
#6:組立場所, 7:待機場所,
#[ 8:L, 8:B, 8:W],[[9:取付],[9:溶接],[9:艤装]],[[10:CL],[10:CA],[10:TB]],
#11:ID1, 12:ID2, 13:due],\          
data={
"A00100":["S1777_3AS1_A00100_8Lx7B_入荷",10,1,188 ,\
["A00101","A00102","A00103","0","0"],[[0,0],[ 0,0],[ 0,0],[ 0,0],[ 0,0]], \
0, 0,\
[ 8 , 7 , 15 ],[[1 ],[0 ],[0 ]],[[0 ],[0 ],[0 ]],\
2, 1, ],\
"A00101":["S1777_3AS1_A00101_8Lx7B_D7_3AS1",73,18,189 ,\
["A00110","0","0","0","0"],[[0,0],[ 0,0],[ 0,0],[ 0,0],[ 0,0]],\
 D7_3AS1, 0,\
[ 8 , 7 , 15 ],[[56 ],[0 ],[0 ]],[[0 ],[0 ],[0 ]],\
2, 1, ],\
 :
}

#プログラムI-2.1b
#=====データ修正
data["A02201"][8][0]=14
data["A06950"][2]=1
data["A07050"][2]=1
data["A07540"][2]=1 
data["A08050"][2]=1
data["A08520"][2]=1 
data["A09030"][2]=1
data["A02610"][5][1][1]=data["A02610"][5][0][1]
data["A03010"][5][1][1]=data["A03010"][5][0][1]
data["A03210"][5][1][1]=data["A03210"][5][0][1]
data["A04310"][5][1][1]=data["A04310"][5][0][1]
data["A09210"][5][0][1]=-3.5
data["A09210"][5][1][1]=data["A09210"][5][0][1]
data["A09322"][4][0],data["A09322"][4][1]=\
data["A09322"][4][1],data["A09322"][4][0]
data0=data
#sys.exit() 

●研究対象とする造船所の全景と物流を次に示します。敷地の制約から、艤装と待機は台船上で行われます。したがって、上構部を除いて、総組は行われていません。

●組立場所のレイアウトを次に示します。

●艤装・待機場所のレイアウトを次に示します。

●カレンダーとしては、いわゆる工場カレンダーのほかに、ブラスト専用のカレンダーが使われています。当研究では、工場カレンダーに基づいて計画単位を半日としています。

工場カレンダー	ブラストカレンダー

●主に日程計画をまとめた総合日程表が作成されており、これを基にして定盤計画と配員計画が実施されています。

●また配員可能数については次を参考にしています。

プログラム：Part I-1

#プログラムI-1.0
from optseq import *
import math
#=====
print("")
print(' ************************ 艤装工程計画 *************************')
#=====
usuki1=Model()
CS=["CS","CC","SS","SC"]

●プログラムI-1.1

#プログラムI-1.1a
#=======================================================================
#リソースの定義（場所） 
#=======================================================================
#作業場所の名前と矩形サイズ [name,L=L2-L1+1,B=B2-B1+1]
R={
 1:["R[K1]",   42,  13  ],\
 2:["R[K2]",   56,  13  ],\
 3:["R[K3]",   47,  15  ],\
 4:["R[K4S4]", 33,  37+2],\
 5:["R[K5]",   45,  23+2],\
 6:["R[K6S2]", 61,  26],\
 7:["R[AB]",   19,  16  ],\
 8:["R[BB]",   19,  16  ],\
 9:["R[CB]",   19,  16  ],\
10:["R[WB]",   14,  13  ],\
11:["R[D35]",  33+2,18+2],\
12:["R[D5]",   54+2,16+2],\
13:["R[D7]",   72,  22+2],\
14:["R[RD1]",  40+2,14+2],\
15:["R[RD3]",  35+2,17+2],\
16:["R[RD4]",  40+2,15+2],\
17:["R[RD5]",  30+2,12+2],\
18:["R[RD6]",  30+2,12+2],\
19:["R[RD7]",  35+2,15+2],\
20:["R[RD8]",  40+2,15+2],\
}

#プログラムI-1.1b
#-----各場所の長さ方向線分[L1,L2]・幅方向線分[B1,B2]
K1=  [ 1,0, R[1][1],0, R[1][2]]
K2=  [ 2,0, R[2][1],0, R[2][2]]
K3=  [ 3,0, R[3][1],0, R[3][2]]
K4S4=[ 4,0, R[4][1],0, R[4][2]]
K5=  [ 5,0, R[5][1],0, R[5][2]]
K6S2=[ 6,0, R[6][1],0, R[6][2]]
AB = [ 7,0, R[7][1],0, R[7][2]]
BB = [ 8,0, R[8][1],0, R[8][2]]
CB = [ 9,0, R[9][1],0, R[9][2]]
WB = [10,0,R[10][1],0,R[10][2]]
D35= [11,0,R[11][1],0,R[11][2]]
D5=  [12,0,R[12][1],0,R[12][2]]
D7=  [13,0,R[13][1],0,R[13][2]]
RD1= [14,0,R[14][1],0,R[14][2]]
RD3= [15,0,R[15][1],0,R[15][2]]
RD4= [16,0,R[16][1],0,R[16][2]]
RD5= [17,0,R[17][1],0,R[17][2]]
RD6= [18,0,R[18][1],0,R[18][2]]
RD7= [19,0,R[19][1],0,R[19][2]]
RD8= [20,0,R[20][1],0,R[20][2]]

#プログラムI-1.1c
#-----各区画の長さ方向線分[L1,L2]・幅方向線分[B1,B2]
K12= [ 1, 14, 28,         0, R[1][2] ]
K13= [ 1, 28, 42,         0, R[1][2] ]
K14= [ 1, 42, R[1][1],    0, R[1][2] ]
K21= [ 2,  0, 14,         0, R[2][2] ]
K22= [ 2, 14, 28,         0, R[2][2] ]
K23= [ 2, 28, 42,         0, R[2][2] ]
K24= [ 2, 42, R[2][1],    0, R[2][2] ]
K312=[ 3,  0, 24,         0, R[3][2] ]
K33 =[ 3, 24, 47,         0, R[3][2] ]
S4  =[ 4, 11, R[4][1],    0, 12      ]
K411=[ 4,  0, 16,        12, 26      ]
K412=[ 4, 16, R[4][1],   12, 26      ]
K421=[ 4,  0, 16,        26, R[4][2] ]
K422=[ 4, 16, R[4][1],   26, R[4][2] ]
K5U= [ 5,  7, R[5][1],    0, 12      ] 
K5R= [ 5,  0, R[5][1],   10+2, R[5][2] ] 
S2W= [ 6,  1, 22,         1, 13+2    ]
S2E= [ 6, 36, R[6][1],    0, 13+2    ]
K6=  [ 6,  0, R[6][1],   13, R[6][2] ]
D7N =[13,  0, 31,         0, 24      ]
D7S =[13, 38, 63,         0, 24      ]
D7EN=[13,  0, 31,         0, 12      ]
D7ES=[13, 38, 63,         0, 12      ]
D7WN=[13,  0, 31,        13, 24      ]
D7WS=[13, 38, 63,        13, 24      ] 

#プログラムI-1.1d
#-----データセットにおける割付場所
K1K2    =[K12,K13,K14,K21,K22,K23,K24]; 
K1K2_K5U=[K12,K13,K14,K21,K22,K23,K24,K5U]
K3=[K312,K33]; 
K4=[K411,K412,K421,K422,S4,WB];
K5U=[K5U];
K5R=[K5R];
K6=[K6,S2W,S2E]
#AB=[AB]; BB=[BB]; CB=[CB];
ABC=[AB,BB,CB]; 
AB=ABC; BB=ABC; CB=ABC;
RD =[RD1,RD3,RD4,RD5,RD6,RD7,RD8,D35,D5,D7N,D7S,S4] 
RD2=[RD1,RD3,RD4,RD5,RD6,RD7,RD8,D35,D5,D7N,D7S,S4] 

#プログラムI-1.1e
#-----ブロック特有の割付場所の起点・終点 [i,x1,y1,x2,y2]
K1_1FS1=[[ 1,R[1][1]-14, R[1][1],   0, R[1][2]]]
K1_1S6P=[[ 1,R[1][1]-11, R[1][1],   0, R[1][2]]]
K1_1S6S=[[ 1,R[1][1]-22, R[1][1]-11,0, R[1][2]]]
K2_8S6P=[[ 2,R[2][1]-22, R[2][1]-11,0, R[2][2]]]
K2_8S6S=[[ 2,R[2][1]-11, R[2][1],   0, R[2][2]]]
K4_1AS1=[S4]
D5_AH=  [[12,1,25,0,12]]
D7_3AS1=[D7ES]
D7_1AS1P=[D7ES]
D7_1AS1S=[D7ES]
D7_1AS5P=[D7WN]
D7_1AS5S=[D7WN]
D7_2S6P=[D7ES]
D7_2S6S=[D7ES]
D7_AH1P=[D7N]
D7_AH1S=[D7N]
D7_AH2=[D7N]
D7_3AS12=[D7ES]
D7_1AS1P2=[D7ES]
D7_1AS1S2=[D7ES]
D7_1AS5P2=[D7WN]
D7_1AS5S2=[D7WN]
県岸=0

#プログラムI-1.1f
resPL={}
for id in R:
  resPL[id]=0
  if not id in [4,5,6,13]:
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):        
        resPL[id,i,j] =usuki1.addResource(R[id][0]+"[{0:02d}_{1:02d}]"\
                                         .format(i,j), capacity={(0,"inf"):1}) 
  if id==4:
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in  range(0,11) and j in range(0,12):
          resPL[id,i,j] =usuki1.addResource(R[id][0]+"[{0:02d}_{1:02d}]"\
                                           .format(i,j), capacity={(0,"inf"):0})
        else:     
          resPL[id,i,j] =usuki1.addResource(R[id][0]+"[{0:02d}_{1:02d}]"\
                                           .format(i,j), capacity={(0,"inf"):1})
  if id==5:
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in  range(0,7) and j in range(0,10):
          resPL[id,i,j] =usuki1.addResource(R[id][0]+"[{0:02d}_{1:02d}]"\
                                           .format(i,j), capacity={(0,"inf"):0})
        else:     
          resPL[id,i,j] =usuki1.addResource(R[id][0]+"[{0:02d}_{1:02d}]"\
                                           .format(i,j), capacity={(0,"inf"):1})
  if id==6:       
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in [0,22,23,24,25,26,27,28,29,30,31,32,33,34,35] and j in range(0,13):
          resPL[id,i,j] =usuki1.addResource(R[id][0]+"[{0:02d}_{1:02d}]"\
                                           .format(i,j), capacity={(0,"inf"):0})
        elif i in range(1,17) and j in range(0,1):
          resPL[id,i,j] =usuki1.addResource(R[id][0]+"[{0:02d}_{1:02d}]"\
                                           .format(i,j), capacity={(0,"inf"):0})            
        else:     
          resPL[id,i,j] =usuki1.addResource(R[id][0]+"[{0:02d}_{1:02d}]"\
                                           .format(i,j), capacity={(0,"inf"):1})
  if id==13:       
    for i in range(0,R[id][1]):
      for j in range(0,R[id][2]):  
        if i in [31,32,33,34,35,36,37,63,64,65,66,67,68,69,70,71] and j in range(0,24):
          resPL[id,i,j] =usuki1.addResource(R[id][0]+"[{0:02d}_{1:02d}]"\
                                           .format(i,j), capacity={(0,"inf"):0})         
        else:     
          resPL[id,i,j] =usuki1.addResource(R[id][0]+"[{0:02d}_{1:02d}]"\
                                           .format(i,j), capacity={(0,"inf"):1})

●プログラムI-1.2

#プログラムI-1.2    
#=======================================================================
#リソースの定義（作業員）
#=======================================================================
eps1=1.25
resWH={}    #取付28人
r=math.floor(28*4*eps1);
resWH=usuki1.addResource("R[WH]", capacity={(0,"inf"):r})
#-----
resWW={}    #溶接43人
r=math.floor(43*4*eps1);
resWW=usuki1.addResource("R[WW]", capacity={(0,"inf"):r})
#-----
eps2=1.25
resWF1={}   #鉄艤5人
r=math.floor(5*4*eps2);
resWF1=usuki1.addResource("R[WF1]", capacity={(0,"inf"):r})
#-----
resWF2={}   #甲配7人
r=math.floor(7*4*eps2);
resWF2=usuki1.addResource("R[WF2]", capacity={(0,"inf"):r})
#-----
resWF3={}   #機配5人⇒10人
r=math.floor(10*4*eps2);
resWF3=usuki1.addResource("R[WF3]", capacity={(0,"inf"):r})
#-----
resWP={}    #塗装内4人⇒8人
r=math.floor(8*4*eps2);
resWP=usuki1.addResource("R[WP]", capacity={(0,"inf"):r})
#-----
resWP2={}   #塗装外8人
r=math.floor(8*4*eps2);
resWP2=usuki1.addResource("R[WP2]", capacity={(0,"inf"):r})

例題(ex11,12,13,14,15,16,17,22,23,24)

●例題ex11：実務的な機械スケジューリング問題

#ex11

●例題ex12：状態変数

#ex12

●例題ex13：順序依存の段取り時間

#ex13

●例題ex14：貯蔵資源の表現方法

#ex14

●例題ex15：ジョブショップスケジューリング問題のベンチマーク問題例

#ex15

●例題ex16：後ろ詰めの表現方法

#ex16

●例題ex17：リリース時刻

#ex17

●例題ex22：Excel出力用の例題

#ex22

●例題23：モードに依存した段取り時間

#ex23

●例題ex24：納期遅れに対する2次のペナルティ

#ex24