fprintf('(1)系の固有値の計算(Matlabでのコマンド)\n')

	fprintf('パラメータの入力')
	M=0.628;
	e=4.01;
	G=3.18;
	m=0.004735;
	J=0.0001155;
	l=0.2465;
	f=0.0001003;
	g=9.8;
	a11=0;
	a12=0;
	a13=1;
	a14=0;
	a21=0;
	a22=0;
	a23=0;
	a24=1;
	a31=0;
	a41=0;
	S=1/((M+m)*J+M*m*l*l);
	a32=(-1)*m*m*l*l*g*S;
	a33=(-1)*e*(J+m*l*l)*S;
	a34=f*m*l*S;
	a42=m*g*l*(m+M)*S;
	a43=m*l*e*S;
	a44=(-1)*(M+m)*f*S;
	fprintf('パラメータの入力完了\n')
	fprintf('Vector a の各成分の計算\n')
	 a=[a11 a12 a13 a14
	a21 a22 a23 a24
	a31 a32 a33 a34 
	a41 a42 a43 a44]
	fprintf('Vector aの固有値の計算\n')
	fprintf('eig(a)\n')
	eig(a)

	%固有値をMatlabに読みこませておく。

	s1=0;s2=5.2094;s3=-5.3711;s4=-6.4600;

fprintf('(2)極配置法によるフィードバックゲインの決定\n')

	fprintf('Vector b の各成分を入力\n')

	b1=0;
	b2=0;
	b3=G*(J+m*l*l)*S;
	b4=(-1)*m*l*G*S;
	b=[b1
	b2
	b3
	b4]

	fprintf('状態方程式の各係数の入力　dx/dt=ax+bu y=cx+du\n')
	
	c=[1 0 0 0
	0 1 0 0
	0 0 1 0
	0 0 0 1]
	
	d=[0
	0
	0
	0]

	%虚数iを定義する。

	i=sqrt(-1)

	fprintf('極の配置を指定する。システムが安定になるように極を選ぶ。\n')
	fprintf('極配置の例1\n')
	p11=-1+5i;
	p12=-1-5i;
	p13=-23;
	p14=-5;

	p1=[p11 p12 p13 p14]

	fprintf('極配置の例2\n')
	p21=-3+10i;
	p22=-3-10i;
	p23=-5+2i;
	p24=-5-2i;

	p2=[p21 p22 p23 p24]
	
	fprintf('フィードバックゲインKを計算する。\n')
	fprintf('例1のフィードバックゲイン\n')
	fprintf('K1= place(a,b,p1)\n')
	K1= place(a,b,p1)
	fprintf('例2のフィードバックゲイン\n')
	fprintf('K2= place(a,b,p2)\n')
	K2= place(a,b,p2)
	fprintf('レギュレータの極L1,L2,L3,L4を計算する。\n')
	fprintf('例1の極\n')
	eig(a-b*K1)
	fprintf('例2の極\n')
	eig(a-b*K2)

	fprintf('シミュレーション結果を出力する。\n')
	
	fprintf('倒立振り子系の初期値の入力\n')
	x0=[0
	0.5
	1
	0]

	fprintf('サンプリングタイムの設定\n')
	fprintf('t=0.0:0.05:5.0\n') 
	t=0.0:0.05:5.0;

	fprintf('例1のシミュレーション\n')
	fprintf('Y1=initial(a-b*K1,b,c,d,x0,t)\n')
	Y1=initial(a-b*K1,b,c,d,x0,t);
	plot(t,Y1),grid;
	xlabel('time[sec]'), ylabel('θ[rad] or Z[m] ')
	pause
	clc

	fprintf('例2のシミュレーション\n')
	fprintf('Y2=initial(a-b*K2,b,c,d,x0,t)\n')
	Y2=initial(a-b*K2,b,c,d,x0,t);
	plot(t,Y2),grid;
	xlabel('time[sec]'), ylabel('θ[rad] or Z[m] ')
	pause
	clc

fprintf('(3)最適レギュレーターの構成\n')

	fprintf('重み行列rの入力\n')
	fprintf('r=1\n')
	r=1;

	fprintf('レギュレーターの設計、フィードバックゲインkの決定\n')

	q1=[1 0 0 0
	0 1 0 0
	0 0 1 0
	0 0 0 1];

	fprintf('k1=lqr(a,b,q1,r)\n')
	k1=lqr(a,b,q1,r)

	fprintf('ケース1についてのシミュレーション\n')
	fprintf('y1=initial(a-b*k1,b,c,d,x0,t)\n')
	y1=initial(a-b*k1,b,c,d,x0,t);
	plot(t,y1),grid;
	xlabel('time[sec]'), ylabel('θ[rad] or Z[m] ')
	pause
	clc
	
	q2=[100 0 0 0
	0 1 0 0
	0 0 1 0
	0 0 0 1];
	
	fprintf('k2=lqr(a,b,q2,r)\n')
	k2=lqr(a,b,q2,r)

	fprintf('ケース2についてのシミュレーション\n')
	fprintf('y2=initial(a-b*k2,b,c,d,x0,t)\n')
	y2=initial(a-b*k2,b,c,d,x0,t);
	plot(t,y2),grid;
	xlabel('time[sec]'), ylabel('θ[rad] or Z[m] ')
	pause
	clc
	
	q3=[1 0 0 0
	0 100 0 0
	0 0 1 0
	0 0 0 1];
	
	fprintf('k3=lqr(a,b,q3,r)\n')
	k3=lqr(a,b,q3,r)

	fprintf('ケース3についてのシミュレーション\n')
	fprintf('y3=initial(a-b*k3,b,c,d,x0,t)\n')
	y3=initial(a-b*k3,b,c,d,x0,t);
	plot(t,y3),grid;
	xlabel('time[sec]'), ylabel('θ[rad] or Z[m] ')
	pause
	clc
	
	q4=[1 0 0 0
	0 1 0 0 
	0 0 100 0
	0 0 0 1];
	
	fprintf('k4=lqr(a,b,q4,r)\n')
	k4=lqr(a,b,q4,r)

	fprintf('ケース4についてのシミュレーション\n')
	fprintf('y4=initial(a-b*k4,b,c,d,x0,t)\n')
	y4=initial(a-b*k4,b,c,d,x0,t);
	plot(t,y4),grid;
	xlabel('time[sec]'), ylabel('θ[rad] or Z[m] ')
	pause
	clc
	
	q5=[1 0 0 0
	0 1 0 0 
	0 0 1 0
	0 0 0 100];
	
	fprintf('k5=lqr(a,b,q5,r)\n')
	k5=lqr(a,b,q5,r)

	fprintf('ケース5についてのシミュレーション\n')
	fprintf('y5=initial(a-b*k5,b,c,d,x0,t)\n')
	y5=initial(a-b*k5,b,c,d,x0,t);
	plot(t,y5),grid;
	xlabel('time[sec]'), ylabel('θ[rad] or Z[m] ')
	pause
	clc