動的システムのシミュレーション

Ⅶ．動的システムのシミュレーション

システムのモデル

　　問題によっては，システムのモデルを作成し，そのモデルを操作した結果を見てみない限り，解が得られない場合も多くあります．システムのモデルとしては，以下に示すように，多くの方法が存在しますが，ここでは，動的システムを表す方法の一つである微分方程式モデルと待ち行列モデルについて説明します．

微分方程式モデル

　　上図に示した倒立振り子に対する微分方程式を解析的に解くことは不可能です（線形化すれば可能）．そのような場合，コンピュータを使用して数値解を求める（コンピュータ上でシステムをシミュレーションする）必要があります．その際，よく使用されるのが，ルンゲ・クッタ法です．下に示すプログラムは，(1) 式を解くためのものであり，各パラメータは以下のようになっています．

m : 10.0
r : 0.5
k : m
u = -k₁θ　　ただし，k₁ = 2 m r g

　　(1) 式が，1 階の連立微分方程式

に変換されて解かれていることに注意してください．

001	/****************************/
002	/* 倒立振り子の微分方程式   */
003	/*      coded by Y.Suganuma */
004	/****************************/
005	#include <stdio.h>
006	#include <math.h>
007	
008	double a[2][2], b, k1;
009	
010	/****************************/
011	/* 倒立振り子の運動方程式   */
012	/*      time : 時間         */
013	/*      x    : 位置と速度   */
014	/*      dx   : xの微分      */
015	/****************************/
016	void link(double time, double *x, double *dx)
017	{
018		dx[0] = a[0][0] * x[0] + a[0][1] * x[1];
019		dx[1] = a[1][0] * sin(x[0]) + a[1][1] * x[1] + b * k1 * x[0];
020	}
021	
022	/*******************************************/
023	/* ルンゲ・クッタ法  dx/dt=f(t,x)          */
024	/*      time : 現在の時間                  */
025	/*      h : 時間刻み幅                     */
026	/*      x : 現在の状態                     */
027	/*      dx : 微係数（f(t,x):snxで計算する）*/
028	/*      g : 作業域（g[4][n]）              */
029	/*      n : 微分方程式の次数               */
030	/*      snx : 微係数を計算する関数の名前   */
031	/*      return : time+h                    */
032	/*******************************************/
033	double rkg(double time, double h, double *x, double *dx, double **g,
034	           int n, void (*sub)(double, double *, double *))
035	{
036		double h2 = 0.5 * h;
037	
038		(*sub)(time, x, dx);
039		for (int i1 = 0; i1 < n; i1++)
040			g[0][i1] = h * dx[i1];
041	
042		time += h2;
043		for (int i1 = 0; i1 < n; i1++)
044			g[1][i1] = x[i1] + 0.5 * g[0][i1];
045		(*sub)(time, g[1], dx);
046		for (int i1 = 0; i1 < n; i1++)
047			g[1][i1] = h * dx[i1];
048	
049		for (int i1 = 0; i1 < n; i1++)
050			g[2][i1] = x[i1] + 0.5 * g[1][i1];
051		(*sub)(time, g[2], dx);
052		for (int i1 = 0; i1 < n; i1++)
053			g[2][i1] = h * dx[i1];
054	
055		time += h2;
056		for (int i1 = 0; i1 < n; i1++)
057			g[3][i1] = x[i1] + g[2][i1];
058		(*sub)(time, g[3], dx);
059		for (int i1 = 0; i1 < n; i1++)
060			g[3][i1] = h * dx[i1];
061	
062		for (int i1 = 0; i1 < n; i1++)
063			x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0;
064	
065		return time;
066	}
067	
068	int main(void)
069	{
070				// 領域の確保
071		double **wk = new double * [4];
072		for (int i1 = 0; i1 < 4; i1++)
073			wk[i1] = new double [2];
074				// 微分方程式の係数の設定
075		double x[2], dx[2];
076		double g   = 9.8;
077		double pi  = 4.0 * atan(1.0);
078		double ut  = pi / 180.0;
079		double uti = 1.0 / ut;
080		double m   = 10.0;
081		double l   = 1.0;
082		double r   = 0.5 * l;
083		double I   = l * l / 12.0;
084		double k   = m;
085		k1         = 2.0 * m * r * g;
086	
087		a[0][0] = 0.0;
088		a[0][1] = 1.0;
089		a[1][0] = m * r * g / (I + m * r * r);
090		a[1][1] = -k / (I + m * r * r);
091		b       = -1.0 / (I + m * r * r);
092				// 初期設定
093		x[0] = 10.0;
094		x[1] = 0.0;
095				// ラジアン単位に変換
096		x[0] *= ut;
097		x[1] *= ut;
098				// 微分方程式を解く
099		double h     = 0.01;
100		double time  = 0.0;
101		double end_t = 2.999;
102		printf("%f %f\n", time, x[0]*uti);
103		while (time < end_t) {
104			time = rkg(time, h, x, dx, wk, 2, link);
105			printf("%f %f\n", time, x[0]*uti);
106		}
107				// 領域の解法
108		for (int i1 = 0; i1 < 4; i1++)
109			delete [] wk[i1];
110		delete [] wk;
111	
112		return 0;
113	}

（出力）

0.000000 10.000000
0.010000 9.990591
0.020000 9.962861
0.030000 9.917557
0.040000 9.855433
0.050000 9.777243
0.060000 9.683747
0.070000 9.575700
0.080000 9.453861
0.090000 9.318981
0.100000 9.171812
0.110000 9.013095
0.120000 8.843570
0.130000 8.663965
0.140000 8.475001
0.150000 8.277390
0.160000 8.071830
0.170000 7.859010
0.180000 7.639606
0.190000 7.414279
0.200000 7.183677
0.210000 6.948432
0.220000 6.709162
0.230000 6.466468
0.240000 6.220935
0.250000 5.973129
0.260000 5.723601
0.270000 5.472883
0.280000 5.221488
0.290000 4.969912
0.300000 4.718631
0.310000 4.468104
0.320000 4.218769
0.330000 3.971047
0.340000 3.725337
0.350000 3.482022
0.360000 3.241465
0.370000 3.004008
0.380000 2.769977
0.390000 2.539677
0.400000 2.313395
0.410000 2.091401
0.420000 1.873945
0.430000 1.661259
0.440000 1.453558
0.450000 1.251039
0.460000 1.053882
0.470000 0.862251
0.480000 0.676292
0.490000 0.496135
0.500000 0.321897
0.510000 0.153675
0.520000 -0.008445
0.530000 -0.164394
0.540000 -0.314116
0.550000 -0.457568
0.560000 -0.594722
0.570000 -0.725562
0.580000 -0.850085
0.590000 -0.968297
0.600000 -1.080220
0.610000 -1.185884
0.620000 -1.285329
0.630000 -1.378608
0.640000 -1.465779
0.650000 -1.546914
0.660000 -1.622090
0.670000 -1.691394
0.680000 -1.754919
0.690000 -1.812768
0.700000 -1.865047
0.710000 -1.911871
0.720000 -1.953360
0.730000 -1.989639
0.740000 -2.020839
0.750000 -2.047094
0.760000 -2.068543
0.770000 -2.085329
0.780000 -2.097597
0.790000 -2.105497
0.800000 -2.109179
0.810000 -2.108797
0.820000 -2.104506
0.830000 -2.096462
0.840000 -2.084822
0.850000 -2.069746
0.860000 -2.051391
0.870000 -2.029916
0.880000 -2.005479
0.890000 -1.978239
0.900000 -1.948352
0.910000 -1.915975
0.920000 -1.881263
0.930000 -1.844368
0.940000 -1.805443
0.950000 -1.764636
0.960000 -1.722096
0.970000 -1.677967
0.980000 -1.632392
0.990000 -1.585510
1.000000 -1.537459
1.010000 -1.488372
1.020000 -1.438379
1.030000 -1.387609
1.040000 -1.336185
1.050000 -1.284228
1.060000 -1.231854
1.070000 -1.179177
1.080000 -1.126307
1.090000 -1.073349
1.100000 -1.020405
1.110000 -0.967572
1.120000 -0.914946
1.130000 -0.862616
1.140000 -0.810668
1.150000 -0.759185
1.160000 -0.708244
1.170000 -0.657920
1.180000 -0.608283
1.190000 -0.559399
1.200000 -0.511330
1.210000 -0.464135
1.220000 -0.417869
1.230000 -0.372583
1.240000 -0.328322
1.250000 -0.285132
1.260000 -0.243051
1.270000 -0.202116
1.280000 -0.162359
1.290000 -0.123809
1.300000 -0.086493
1.310000 -0.050433
1.320000 -0.015649
1.330000 0.017844
1.340000 0.050032
1.350000 0.080905
1.360000 0.110455
1.370000 0.138677
1.380000 0.165570
1.390000 0.191134
1.400000 0.215372
1.410000 0.238289
1.420000 0.259892
1.430000 0.280191
1.440000 0.299199
1.450000 0.316928
1.460000 0.333394
1.470000 0.348616
1.480000 0.362610
1.490000 0.375399
1.500000 0.387004
1.510000 0.397448
1.520000 0.406756
1.530000 0.414954
1.540000 0.422069
1.550000 0.428129
1.560000 0.433162
1.570000 0.437198
1.580000 0.440268
1.590000 0.442402
1.600000 0.443632
1.610000 0.443990
1.620000 0.443509
1.630000 0.442221
1.640000 0.440159
1.650000 0.437356
1.660000 0.433847
1.670000 0.429663
1.680000 0.424839
1.690000 0.419408
1.700000 0.413403
1.710000 0.406857
1.720000 0.399803
1.730000 0.392273
1.740000 0.384299
1.750000 0.375913
1.760000 0.367146
1.770000 0.358029
1.780000 0.348592
1.790000 0.338865
1.800000 0.328876
1.810000 0.318656
1.820000 0.308230
1.830000 0.297627
1.840000 0.286874
1.850000 0.275994
1.860000 0.265015
1.870000 0.253959
1.880000 0.242851
1.890000 0.231712
1.900000 0.220566
1.910000 0.209432
1.920000 0.198332
1.930000 0.187284
1.940000 0.176307
1.950000 0.165419
1.960000 0.154637
1.970000 0.143977
1.980000 0.133453
1.990000 0.123081
2.000000 0.112873
2.010000 0.102844
2.020000 0.093003
2.030000 0.083364
2.040000 0.073935
2.050000 0.064727
2.060000 0.055748
2.070000 0.047006
2.080000 0.038509
2.090000 0.030263
2.100000 0.022273
2.110000 0.014546
2.120000 0.007085
2.130000 -0.000106
2.140000 -0.007024
2.150000 -0.013666
2.160000 -0.020031
2.170000 -0.026116
2.180000 -0.031922
2.190000 -0.037448
2.200000 -0.042695
2.210000 -0.047663
2.220000 -0.052353
2.230000 -0.056768
2.240000 -0.060910
2.250000 -0.064782
2.260000 -0.068385
2.270000 -0.071725
2.280000 -0.074805
2.290000 -0.077628
2.300000 -0.080200
2.310000 -0.082525
2.320000 -0.084608
2.330000 -0.086455
2.340000 -0.088071
2.350000 -0.089461
2.360000 -0.090632
2.370000 -0.091591
2.380000 -0.092342
2.390000 -0.092893
2.400000 -0.093251
2.410000 -0.093421
2.420000 -0.093411
2.430000 -0.093228
2.440000 -0.092879
2.450000 -0.092369
2.460000 -0.091708
2.470000 -0.090900
2.480000 -0.089955
2.490000 -0.088878
2.500000 -0.087676
2.510000 -0.086357
2.520000 -0.084928
2.530000 -0.083395
2.540000 -0.081764
2.550000 -0.080044
2.560000 -0.078240
2.570000 -0.076359
2.580000 -0.074408
2.590000 -0.072392
2.600000 -0.070318
2.610000 -0.068192
2.620000 -0.066019
2.630000 -0.063807
2.640000 -0.061559
2.650000 -0.059283
2.660000 -0.056983
2.670000 -0.054664
2.680000 -0.052331
2.690000 -0.049990
2.700000 -0.047644
2.710000 -0.045299
2.720000 -0.042959
2.730000 -0.040628
2.740000 -0.038310
2.750000 -0.036008
2.760000 -0.033727
2.770000 -0.031469
2.780000 -0.029239
2.790000 -0.027039
2.800000 -0.024873
2.810000 -0.022742
2.820000 -0.020650
2.830000 -0.018599
2.840000 -0.016591
2.850000 -0.014629
2.860000 -0.012714
2.870000 -0.010848
2.880000 -0.009033
2.890000 -0.007269
2.900000 -0.005560
2.910000 -0.003904
2.920000 -0.002305
2.930000 -0.000762
2.940000 0.000724
2.950000 0.002153
2.960000 0.003523
2.970000 0.004834
2.980000 0.006087
2.990000 0.007280
3.000000 0.008415

待ち行列モデル

　　ここでは，上図に示すような待ち行列システムについて検討します．このシステムは，待ち行列システムとして最も簡単な例です．客が到着すると，窓口が空いていればサービスを受け，空いていなければ唯一存在する待ち行列に並ぶというシステムです．到着分布が母数 λ の指数分布，サービス分布が母数 μ の指数分布に従い，窓口の数が s であり，かつ，システムの容量が無限大であるシステムを対象として考えてみます．なお，指数分布は，ポアソン分布と強い関係があります．例えば，平均到着間隔が平均 1 / λ の指数分布をするとき，単位時間内に到着する到着人数の分布は平均 λ のポアソン分布をします．

　　到着分布は，母数 λ の指数分布ですので，その確率密度関数と分布関数は以下のようになります．

密度関数　f(x) = λe^-λx　　x ≧ 0

　　　　　　　　　 = 0　　　　 x < 0

分布関数　F(x) = 1 - e^-λx　　x ≧ 0　　　x 時間内に客が到着する確率

　　　　　　　　　 = 0　　　　　　 x < 0

このとき，母数 λ は，到着率（単位時間に到着する平均客数）と呼ばれ，また， 1 / λ は，平均到着時間間隔に相当します．

　　また，サービス分布は，母数 μ の指数分布ですので，その確率密度関数と分布関数は以下のようになります．

密度関数　f(x) = μe^-μx　　x ≧ 0

　　　　　　　　　 = 0　　　　 x < 0

分布関数　F(x) = 1 - e^-μx　　x ≧ 0　　　x 時間内にサービスが終了する確率

　　　　　　　　　 = 0　　　　　　 x < 0

このとき，母数 μ は，サービス率（単位時間に窓口が処理する平均客数）と呼ばれ，また， 1 / μ は，平均サービス時間に相当します．

　　微分方程式モデルのように，動的なシステムをシミュレーションする場合は，コンピュータ内部で時間を進めていく必要があります．微分方程式モデルでは，刻み幅という形で，一定時間毎に時間を進めていきました．しかし，待ち行列システムの場合はどうでしょうか．例えば，自動券売機を何台設置すべきかをシミュレーションによって決めたいとします．このとき，時間を進める基本単位である刻み幅をどの程度にすればよいでしょうか．かなり小さくしないと，たくさんの人がほぼ同時に券売機に殺到したような場合に，適切な処理を行うことができません．しかし，誰も来ないときは，システムの状態が変化しないため計算をする必要が無いにもかかわらず，無駄な計算を頻繁に行わなければ成りません．

　　待ち行列システムのように，何か事象（客の到着，サービスの開始・終了等）が発生しない限り状態が変化しないようなシステムを事象駆動システムと呼びます．このようなシステムでは，時間を一定毎に進めず，事象が発生した時点へ次々と時間を進めていく方法がとられます．つまり，次に発生する事象の内，最も早く発生する事象を調べ，その時点まで時間を進め状態の変更等何らかの処理を行うといった方法です．この方法に従い，待ち行列システムをシミュレーションするための基本的な流れを以下に示します．実際のプログラムに関しては，その下に示すプログラム例を参照して下さい．

001	/******************************/
002	/* 簡単な待ち行列問題（M/M/s）*/
003	/*      coded by Y.Suganuma   */
004	/******************************/
005	#include <stdio.h>
006	#include <stdlib.h>
007	#include <ctime>
008	#include <math.h>
009	#include <map>
010	#include <queue>
011	#include "MT.h"
012	using namespace std;
013	
014	/************************/
015	/* クラスCustomerの定義 */
016	/************************/
017	class Customer
018	{
019		public :
020			double time;   // 到着時刻
021			int state;   // 客の状態
022	                     //   =-1 : 待ち行列に入っている
023	                     //   >=0 : サービスを受けている窓口番号
024			/*********************/
025			/* コンストラクタ    */
026			/*      s : 状態     */
027			/*      t : 到着時刻 */
028			/*******************************/
029			Customer::Customer(int s, double t)
030			{
031				time  = t;
032				state = s;
033			}
034	};
035	
036	/**********************/
037	/* クラスQ_baseの定義 */
038	/**********************/
039	class Q_base {
040			int s;   // 窓口の数
041			int asb;   // 全窓口の空き状況
042	                   //    =0 : すべての窓口がふさがっている
043	                   //    =n : ｎ個の窓口が空いている
044			int *sb;   // 各窓口の空き状況
045	                   //    =0 : 窓口は空いている
046	                   //    >0 : サービスを受けている客番号
047			double asb_t;   // すべての窓口がふさがった時刻
048			double c_asb;   // すべての窓口がふさがっている時間の累計
049			double *c_sb;   // 各窓口がふさがっている時間の累計
050			double *st_e;   // 各窓口のサービス終了時刻
051			double *st_s;   // 各窓口がふさがった時刻
052			int m_lq;   // 最大待ち行列長
053			double c_lq_t;   // 待ち行列長にその長さが継続した時間を乗じた値の累計
054			double c_wt;   // 待ち時間の累計
055			double lq_t;   // 現在の待ち行列長になった時刻
056			double m_wt;   // 最大待ち時間
057			double c_sc_t;   // 系内客数にその数が継続した時間を乗じた値の累計
058			double c_sys;   // 滞在時間の累計
059			double m_sys;   // 最大滞在時間
060			double sc_t;   // 現在の系内客数になった時刻
061			int m_sc;   // 最大系内客数
062			double m_a;   // 到着時間間隔の平均値
063			double m_s;   // サービス時間の平均値
064			double at;   // 次の客の到着時刻（負の時は，終了）
065			double p_time;   // 現在時刻
066			double end;   // シミュレーション終了時刻
067			int nc;   // 到着客数カウンタ
068			map<int, Customer> cus;   // 系内にいる客のリスト
069			queue<int> que;   // 待ち行列
070		public:
071			Q_base(int, double, double, double);   // コンストラクタ
072			~Q_base();   // デストラクタ
073			double Next_at();   // 次の到着時刻
074			double Next_sv();   // サービス終了時刻
075			void Control();   // 全体の制御
076			int Next();   // 次の処理の決定
077			int End_o_s();   // 終了処理
078			void Arrive();   // 客の到着処理
079			void Service(int);   // サービス終了
080			void Output();   // 出力
081	};
082	
083	/*****************************************/
084	/* コンストラクタ                        */
085	/*      s_i : 窓口の数                   */
086	/*      m_a_i : 到着時間間隔の平均値     */
087	/*      m_s_i : サービス時間の平均値     */
088	/*      end_i : シミュレーション終了時刻 */
089	/*****************************************/
090	Q_base::Q_base (int s_i, double m_a_i, double m_s_i, double end_i)
091	{
092	/*
093	          設定
094	*/
095		s   = s_i;
096		m_a = m_a_i;
097		m_s = m_s_i;
098		end = end_i;
099	/*
100	          領域の確保
101	*/
102		sb   = new int [s];
103		c_sb = new double [s];
104		st_e = new double [s];
105		st_s = new double [s];
106	
107		for (int i1 = 0; i1 < s; i1++) {
108			sb[i1]   = 0;
109			c_sb[i1] = 0.0;
110		}
111	/*
112	          初期設定
113	*/
114		p_time = 0.0;
115		nc     = 0;
116		asb    = s;
117		m_lq   = 0;
118		m_sc   = 0;
119		c_asb  = 0.0;
120		c_wt   = 0.0;
121		m_wt   = 0.0;
122		c_lq_t = 0.0;
123		lq_t   = 0.0;
124		m_sys  = 0.0;
125		c_sys  = 0.0;
126		c_sc_t = 0.0;
127		sc_t   = 0.0;
128	/*
129	          乱数の初期設定
130	*/
131		init_genrand((unsigned)time(NULL));
132	/*
133	          最初の客の到着時刻の設定
134	*/
135		at = p_time + Next_at();
136	}
137	
138	/****************/
139	/* デストラクタ */
140	/****************/
141	Q_base::~Q_base()
142	{
143		delete [] sb;
144		delete [] c_sb;
145		delete [] st_e;
146		delete [] st_s;
147	}
148	
149	/********************************/
150	/* 次の客の到着までの時間の発生 */
151	/********************************/
152	double Q_base::Next_at()
153	{
154		return -m_a * log(genrand_real3());
155	}
156	
157	/************************/
158	/* サービス時間の発生   */
159	/************************/
160	double Q_base::Next_sv()
161	{
162		return -m_s * log(genrand_real3());
163	}
164	
165	/**************/
166	/* 全体の制御 */
167	/**************/
168	void Q_base::Control()
169	{
170		int sw = 0;
171		while (sw > -2) {
172			sw = Next();   // 次の処理の選択
173			if (sw == -1)
174				sw = End_o_s();   // シミュレーションの終了
175			else {
176				if (sw == 0)
177					Arrive();   // 客の到着処理
178				else
179					Service(sw);   // サービスの終了
180			}
181		}
182	}
183	
184	/**************************************************/
185	/* 次の処理の決定                                 */
186	/*      return : =-1 : シミュレーションの終了     */
187	/*               =0  : 客の到着処理               */
188	/*               =i  : ｉ番目の窓口のサービス終了 */
189	/**************************************************/
190	int Q_base::Next()
191	{
192		int sw = -1;
193		double t  = end;   // シミュレーション終了時刻で初期設定
194						// 次の客の到着時刻
195		if (at >= 0.0 && at < t) {
196			sw = 0;
197			t  = at;
198		}
199						// サービス終了時刻
200		for (int i1 = 0; i1 < s; i1++) {
201			if (sb[i1] > 0 && st_e[i1] <= t) {
202				sw = i1 + 1;
203				t  = st_e[i1];   // 窓口i1のサービス終了時刻
204			}
205		}
206	
207		return sw;
208	}
209	
210	/**********************************/
211	/* 終了処理                       */
212	/*      return : =-1 : 終了前処理 */
213	/*               =-2 : 実際の終了 */
214	/**********************************/
215	int Q_base::End_o_s()
216	{
217		int sw = -2;
218		p_time = end;   // 現在時刻
219		at     = -1.0;   // 次の客の到着時刻
220	
221		for (int i1 = 0; i1 < s; i1++) {
222			if (sb[i1] > 0) {   // サービス中の客はいるか？
223				if (sw == -2) {
224					sw  = -1;
225					end = st_e[i1];   // 窓口i1のサービス終了時刻
226				}
227				else {
228					if (st_e[i1] > end)
229						end = st_e[i1];   // 窓口i1のサービス終了時刻
230				}
231			}
232		}
233	
234		return sw;
235	}
236	
237	/****************/
238	/* 客の到着処理 */
239	/****************/
240	void Q_base::Arrive()
241	{
242	/*
243	          客数の増加と次の客の到着時刻の設定
244	*/
245		nc     += 1;   // 到着客数カウンタ
246		p_time  = at;   // 現在時刻
247		at      = p_time + Next_at();      // 次の客の到着時刻
248		if (at >= end)
249			at = -1.0;
250	/*
251	          系内客数の処理
252	*/
253		c_sc_t += cus.size() * (p_time - sc_t);   // 系内客数にその数が継続した時間を乗じた値の累計
254		sc_t    = p_time;   // 現在の系内客数になった時刻
255		if ((int)cus.size()+1 > m_sc)
256			m_sc = cus.size() + 1;   // 最大系内客数
257	/*
258	          窓口に空きがない場合
259	*/
260		if (asb == 0) {
261			Customer ct_p(-1, p_time);
262			cus.insert(pair<int, Customer>(nc, ct_p));   // 客の登録（系内客数）
263			c_lq_t += que.size() * (p_time - lq_t);   // 待ち行列長にその長さが継続した時間を乗じた値の累計
264			que.push(nc);   // 客の登録（待ち行列）
265			lq_t    = p_time;   // 現在の待ち行列長になった時刻
266			if ((int)que.size() > m_lq)
267				m_lq = que.size();   // 最大待ち行列長
268		}
269	/*
270	          すぐサービスを受けられる場合
271	*/
272		else {
273			int k = -1;
274			for (int i1 = 0; i1 < s && k < 0; i1++) {
275				if (sb[i1] == 0) {
276					Customer ct_p(i1, p_time);
277					cus.insert(pair<int, Customer>(nc, ct_p));   // 客の登録（系内客数）
278					k        = i1;
279					sb[k]    = nc;   // サービスを受けている客番号
280					st_e[k]  = p_time + Next_sv();   // 窓口kのサービス終了時刻
281					asb     -= 1;   // 空いている窓口数
282				}
283			}
284			st_s[k] = p_time;   // 窓口kがふさがった時刻
285			if (asb == 0)
286				asb_t = p_time;   // すべての窓口がふさがった時刻
287		}
288	}
289	
290	/*********************************/
291	/* サービス終了時の処理          */
292	/*      k : サービス終了窓口番号 */
293	/*********************************/
294	void Q_base::Service(int k)
295	{
296	/*
297	          時間の設定
298	*/
299		k      -= 1;
300		p_time  = st_e[k];   // 現在時刻
301		st_e[k] = -1.0;   // サービス終了時間
302	/*
303	          系内客数の処理
304	*/
305		c_sc_t += cus.size() * (p_time - sc_t);   // 系内客数にその数が継続した時間を乗じた値の累計
306		sc_t    = p_time;   // 現在の系内客数になった時刻
307	/*
308	          滞在時間の処理
309	*/
310		map<int, Customer>::iterator it = cus.find(sb[k]);
311		double x1 = p_time - (it->second).time;
312		c_sys += x1;   // 滞在時間の累計
313		if (x1 > m_sys)
314			m_sys = x1;   // 最大滞在時間
315		cus.erase(sb[k]);   // 客の削除（系内客数）
316	/*
317	          窓口使用時間の処理
318	*/
319		asb     += 1;   // 空いている窓口数
320		sb[k]    = 0;   // 窓口kを空き状態にする
321		c_sb[k] += (p_time - st_s[k]);   // 窓口kがふさがっている時間の累計
322	/*
323	          待ち行列がある場合
324	*/
325		if (que.size() > 0) {
326			asb    -= 1;   // 開いている窓口数
327			c_lq_t += que.size() * (p_time - lq_t);   // 待ち行列長にその長さが継続した時間を乗じた値の累計
328			int n = que.front();
329			que.pop();   // 客の削除（待ち行列）
330			lq_t    = p_time;   // 現在の待ち行列長になった時刻
331			it      = cus.find(n);
332			x1      = p_time - (it->second).time;
333			c_wt   += x1;   // 待ち時間の累計
334			if (x1 > m_wt)
335				m_wt = x1;   // 最大待ち時間
336			int k = -1;
337			for (int i1 = 0; i1 < s && k < 0; i1++) {
338				if (sb[i1] == 0) {
339					k        = i1;
340					sb[k]    = n;   // 窓口kの客番号
341					st_e[k]  = p_time + Next_sv();   // 窓口kのサービス終了時刻
342					st_s[k]  = p_time;   // 窓口kがふさがった時刻
343					(it->second).state = k;   // 客の状態変更
344				}
345			}
346		}
347	/*
348	          待ち行列がない場合
349	*/
350		else {
351			if (asb == 1)
352				c_asb += (p_time - asb_t);   // すべての窓口がふさがっている時間の累計
353		}
354	}
355	
356	/************************/
357	/* 結果の出力           */
358	/* （カッコ内は理論値） */
359	/************************/
360	void Q_base::Output()
361	{
362		double rn  = (double)nc;
363		double rs  = (double)s;
364		double ram = 1.0 / m_a;
365		double myu = 1.0 / m_s;
366		double rou = ram / (rs * myu);
367		double p0, pi;
368		if (s == 1) {
369			p0 = 1.0 - rou;
370			pi = rou;
371		}
372		else {
373			p0 = 1.0 / (1.0 + 2.0 * rou + 4.0 * rou * rou / (2.0 * (1.0 - rou)));
374			pi = 4.0 * rou * rou * p0 / (2.0 * (1.0 - rou));
375		}
376		double Lq = pi * rou / (1.0 - rou);
377		double L  = Lq + rs * rou;
378		double Wq = Lq / ram;
379		double W  = Wq + 1.0 / myu;
380		printf("シミュレーション終了時間＝%.3f 客数＝%d ρ＝%.3f ｐ0＝%.3f\n",
381	           p_time, nc, rou, p0);
382		printf("   すべての窓口が塞がっている割合＝%.3f (%.3f)\n",
383	           c_asb/p_time, pi);
384		printf("   各窓口が塞がっている割合\n");
385		for (int i1 = 0; i1 < s; i1++)
386			printf("      %d   %.3f\n", i1+1, c_sb[i1]/p_time);
387		printf("   平均待ち行列長＝%.3f (%.3f)  最大待ち行列長＝%d\n",
388	           c_lq_t/p_time, Lq, m_lq);
389		printf("   平均系内客数  ＝%.3f (%.3f)  最大系内客数  ＝%d\n",
390	           c_sc_t/p_time, L, m_sc);
391		printf("   平均待ち時間  ＝%.3f (%.3f)  最大待ち時間  ＝%.3f\n",
392	           c_wt/rn, Wq, m_wt);
393		printf("   平均滞在時間  ＝%.3f (%.3f)  最大滞在時間  ＝%.3f\n",
394	           c_sys/rn, W, m_sys);
395	}
396	
397	/****************/
398	/* main program */
399	/****************/
400	int main()
401	{
402		int s;
403		double end, m_a, m_s;
404	/*
405	          入力データ
406	*/
407		printf("窓口の数は？ ");
408		scanf("%d", &s);
409	
410		printf("シミュレーション終了時間？ ");
411		scanf("%lf", &end);
412	
413		printf("   到着時間間隔の平均値は？ ");
414		scanf("%lf", &m_a);
415	
416		printf("   サービス時間の平均値は？ ");
417		scanf("%lf", &m_s);
418	/*
419	          システムの定義
420	*/
421		Q_base base(s, m_a, m_s, end);
422	/*
423	          シミュレーションの実行
424	*/
425		base.Control();
426	/*
427	          出力
428	*/
429		base.Output();
430	
431		return 0;
432	}

（出力）（赤字は入力値であり，括弧内に書かれた数値は理論値）

窓口の数は？ 2
シミュレーション終了時間？ 100000
   到着時間間隔の平均値は？ 3
   サービス時間の平均値は？ 4
シミュレーション終了時間＝100000.116 客数＝33458 ρ＝0.667 ｐ0＝0.200
   すべての窓口が塞がっている割合＝0.543 (0.533)
   各窓口が塞がっている割合
      1   0.726
      2   0.619
   平均待ち行列長＝1.107 (1.067)  最大待ち行列長＝19
   平均系内客数  ＝2.452 (2.400)  最大系内客数  ＝21
   平均待ち時間  ＝3.309 (3.200)  最大待ち時間  ＝43.432
   平均滞在時間  ＝7.329 (7.200)  最大滞在時間  ＝64.841

［演習７］様々な待ち行列システムと考えられる具体例を探し，上で示したプログラム（プログラムを修正しても良い）によってシミュレーションを行い，その結果を報告せよ．