改訂日：2014 年 4 月 18 日

Arduino と MATLAB で制御系設計をはじめよう！
「2012 年１１月 30 日 初版第１刷」 改訂表
改訂箇所頁v

目次P006１.４

本文

改 訂 内 容（旧）1.5

本書で作成したソースコードについて

（新）

1.5

本書で作成したソースコードと補足情報について次の文章に書き改めました。
本書執筆時の動作環境を表1.1に示します。MATLAB/Simulink のバージョンについて，RoTH を試す場合はR2012a 以降が必要です。また，Student Version での動作も確認しています。なお，ArduinoIOについては，バージョンアップにともなって古いバージョンは非公開になるようです。2014 年3 月現在の最新版は4.4 ですが，問題無く動作することを確認しています。
表1.1 本書執筆時の作動環境
ホストPC のOS

Windows 7 32bit

MATLAB/Simulink のバージョン

R2012aMATLAB/Simulink のツールボックス Control System Toolbox （必須）

P006

１.５

見出しP006
本文

１.５

ArduinoUno R3 またはMega 2560 R3

Arduino IDE のバージョン

Version 1.01

ArduinoIO のバージョン

Version 4.2（旧）本書で作成したソースコードについて（新）

本書で作成したソースコードと補足情報について次の文章に書き改めました。
本書で作成したMATLAB のm-file やSimulink モデルについては，下記のURL からダウンロードできます。また，最新のMATLAB/Simulink での動作状況といった補足情報も掲載していますので，是非ご参照ください。
http://books.techshare.co.jp/

P006

１.６

本文

次の文章に書き改めました。
本書で製作した実験装置は，街のホビーショップなどで一般的に入手できる部材を使っています。しかしながら，第5 章で製作するモータ角度制御実験装置においてのみ，角度検出のポテンショメータを取り付けるための専用プレートを使っています。アルミ板等を加工して製作できますが，加工が苦手な方のために，実験キットを用意しました。詳しくは下記Web サイトをご覧ください。
http://store.techshare.jp/

P023

２.５.１

３行目

P025

２.５.３

４行目

P025

補足欄

（誤）

本書執筆時点のバージョンは4.4です。

（正）

本書執筆時点のバージョンは4.2です。

（旧）…次のスケッチを開きます。

（新）…次のスケッチを開きます２。

追加

２

ArduinoIO のバージョンによって異なります。2014年3月現在，最新版（4.4）ではC:\Arduino\ArduinoIO\pde\adioes\adioes.pde になります。
P026

P026

P026

下から４行目下から３行目補足欄（旧）…モータシールドに「Adafruit Motor/Stepper/Servo Shield」２が、…

（新）…モータシールドに「Adafruit Motor/Stepper/Servo Shield」３が、…

（旧）…そのためのAFMotorライブラリ３を…

（新）…そのためのAFMotorライブラリ４を…

（旧）２

http://www.adafruit.com/products/81３

https://github.com/adafruit/Adafruit-Motor-Shield-library/zipball/master（新）３

http://www.adafruit.com/products/81

４

https://github.com/adafruit/Adafruit-Motor-Shield-library/zipball/master

P045

2 行目

（誤）

a.digitalRead(0);（正）a.analogRead(0);

P048

P055

補足欄

４.２

７行目

（誤）

現在のバージョン（4.4）では、…

（正）

本書執筆時のバージョン（4.2）では、…

（旧）

また、実態配線図を図4.5に示します。

（新）また，モータとタコメータにはマブチモータ製FA-130RA-2270を使用しました。実体配線図を図4.5に示します。

P056

図 4.3

2SK2231 とモータとの結線に誤りがありました。
（誤）

+5V
IOREF
RESET
3.3V
5V
GND
GND
VinA0A1A2A3A4A5

AREF
GND
131211109Arduino
87654321+4.5V
1.5kΩ
DG

Motor
S

2SK2231+5V

FET amplifier
10kΩ
100kΩ

M
1N4001

0.1μF

M

指令値入力用
ポテンショ
0.1μF

1N4001

0.1μF

Low pass filter

（正）

P057

図４.５

2SK2231 とモータとの結線に誤りがありました。
（誤）

（正）

P090

プログラム 4-4

⇒付属のソースコードをご参照下さい。

P097

補足欄

（旧）

下から３行目

専用プレートを含む製作キットが TechShare Store から販売される予定です。

（新）専用プレートを含む実験キットがTechShare Store から販売されています。(http://store.techshare.jp/)

P106

下から 2 行目

P111 プログラム 5-3
P124

P126

9 行目

１行目

（誤）

ζ=1の不足減衰では、…

（正）

ζ= 0.2の不足減衰では、…

⇒付属のソースコードをご参照下さい。
（誤）これらの式をKP,KIについて…

（正）

これらの式をKP,KDについて…

（旧）omega_n=10, zeta=0.6に…

（新）

一例として、omega_n=10, zeta=0.6に…

P170

プログラム 6-3

⇒付属のソースコードをご参照下さい。

P175

プログラム 6-4

⇒付属のソースコードをご参照下さい。

P176

プログラム 6-4

⇒付属のソースコードをご参照下さい。

P190

プログラム 7-1

⇒付属のソースコードをご参照下さい。

P199

プログラム 7-2

⇒付属のソースコードをご参照下さい。

P200

プログラム 7-1

⇒付属のソースコードをご参照下さい。

4.6 極指定法による PI ゲイン設計 (p.90)
プログラム 4-4 [velo pi mbd.m]1234567891011121314151617181920212223242526272829%% velo_pi_mbd.m%% Initialize & load dataclose allclear allload model_data%% PI design by pole placement
p1 = input(’p1 = ’);p2 = p1; % 重根の場合% p2 = input(’p2 = ’); % p2 も指定する場合
Kp = -((p1+p2)*T + 1)/K;
Ki = p1*p2*T/K;%% Display PI parametersdisp(’>>> PI parameters <<<’)
fprintf(’Kp = %f\n’,Kp);
fprintf(’Ki = %f\n’,Ki);%% Open simulink model
open_system(’velo_pi_mbd_sl’);
open_system(’velo_pi_mbd_sl/Output’);%% Experiment
ts = 1/50;
sim(’velo_pi_mbd_sl’)%% EOF of velo_pi_mbd.m■ 5.3.3

ステップ応答実験 (p.111)
プログラム 5-3 [pos id step.m]

123456789101112
13
1415
16
1718
19202122232425
262728
2930
3132
33
34
35
3637383940
41424344%% pos_id_step.m%% Initialize
close allclear allload sim_param%% Parameters for identification
r= 60;
r_cyc = 8;
Kp_id = 0.12;
Ncyc= 2;
tfinal = r_cyc*Ncyc;%% ID Experiment
open_system(’pos_id_step_sl’)
open_system(’pos_id_step_sl/Scope’)sim(’pos_id_step_sl’)%% Data processingy = yout.signals(1).values(:,2);t = yout.time;
NN = length(y);N = r_cyc/ts;
yy = reshape(y(2:NN),N,(NN-1)/N);
yf = yy(1:N/2,2:end); % 最初のデータを除き＆前進データのみとする% 平均化と正規化処理
ym = mean(yf’)’;y0 = ym(1); yN = ym(end);ym = (ym-y0)/(yN-y0);%% Plot figuret = (0:N/2-1)*ts;
figure(1)subplot(211)
plot(t,yf), gridxlabel(’Time [s]’),ylabel(’Output [deg]’)
subplot(212)plot(t,ym), gridxlabel(’Time [s]’),ylabel(’Output [deg]’)%% EOF of pos_id_step.m

■ 6.4.2

I-PD+FF 制御のための Simulink モデル (p.170)
プログラム 6-3 [bb pid rc.m]

1234
5
678
910
11
1213
14
15
16
1718
1920
21222324
25
262728
2930
3132
33
34
35
3637
38
3940
41
42
434445464748495051525354555657%% bb_pid_rc.m%% Initialize & load dataclose allclear allload sim_param%% サーボ 1 度 あたりのアームの傾き% RC サーボの場合
K_theta_rc = (pi/180)*(1.5/11.5); % [rad/deg]%% ビーム傾き角 [rad] -> ボール位置 [cm] までの 1/s^2 のゲイン
K_b = (3/5*9.8)*100;
Pb = K_b/s^2;
Pbd = c2d(Pb,ts,’zoh’);[numbd,denbd] = tfdata(Pbd,’v’);%% PID パラメータ for Ball 位置制御
omega_n = 1.5;
zeta= 0.6;
alpha= 0.5;
p1 = (-zeta + j*sqrt(1-zeta^2))*omega_n;p2 = (-zeta - j*sqrt(1-zeta^2))*omega_n;p3 = -alpha;kp = (p1*p2 + p2*p3 + p3*p1)/K_b;
ki = -p1*p2*p3/K_b;
kd = -(p1+p2+p3)/K_b;% kff = 0;kff = ki/alpha;disp(’>>> PID parameters for Ball and Beam <<<’)
fprintf(’kp = %f\n’,kp);
fprintf(’ki = %f\n’,ki);
fprintf(’kd = %f\n’,kd);
fprintf(’kff = %f\n’,kff);%% Reference for ball position
r1 = 20-5;
r2 = 20+5;%% LPF カットオフ周波数
wf = 2*pi*2;%% LPF for ball position sensor
Fc = wf^2/(s^2+2*0.7*wf*s+wf^2);
Fd = c2d(Fc,ts,’tustin’);[numlpf,denlpf] = tfdata(Fd,’v’);%% Open simulink model
open_system(’bb_pid_rc_sl’);
open_system(’bb_pid_rc_sl/Ref angle’)
open_system(’bb_pid_rc_sl/Ball position and reference’)%% EOF of bb_pid_rc.m■ 6.5.3

実験してみよう (p.175–176)
プログラム 6-4 [bb pid hg.m]12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273%% bb_pid_hg.m%% Initialize & load dataclose allclear allload sim_param%% Set identified parameters
K = 300;
T = 0.28;%% Discrete-time plant model
P = K/(T*s^2 + s);
Pd = c2d(P,ts,’zoh’);[numpd,denpd] = tfdata(Pd,’v’);%% PD パラメータ for 自作サーボ
omega_n_hg = 10;
zeta_hg= 0.65;p1 = (-zeta_hg + j*sqrt(1-zeta_hg^2))*omega_n_hg;
p2 = (-zeta_hg - j*sqrt(1-zeta_hg^2))*omega_n_hg;%% Set PD parameters
Kp = p1*p2*T/K;
Kd = -((p1+p2)*T + 1)/K;
Ki = 0;%% Display PID parametersdisp(’>>> PID parameters for HG Servo <<<’)
fprintf(’Kp = %f\n’,Kp);
fprintf(’Ki = %f\n’,Ki);
fprintf(’Kd = %f\n’,Kd);%% サーボ 1 度 あたりのアームの傾き% 自作サーボの場合
K_theta_hg = (pi/180)*(2.1/15); % [rad/deg]%% ビーム傾き角 [rad] -> ボール位置 [cm] までの 1/s^2 のゲイン
K_b = (3/5*9.8)*100;
Pb = K_b/s^2;
Pbd = c2d(Pb,ts,’zoh’);[numbd,denbd] = tfdata(Pbd,’v’);%% PID パラメータ for Ball 位置制御
omega_n = 1.5;
zeta
= 0.6;
alpha= 0.5;
p1 = (-zeta + j*sqrt(1-zeta^2))*omega_n;p2 = (-zeta - j*sqrt(1-zeta^2))*omega_n;p3 = -alpha;kp = (p1*p2 + p2*p3 + p3*p1)/K_b;
ki = -p1*p2*p3/K_b;
kd = -(p1+p2+p3)/K_b;% kff = 0;kff = ki/alpha;disp(’>>> PID parameters for Ball and Beam <<<’)
fprintf(’kp = %f\n’,kp);
fprintf(’ki = %f\n’,ki);
fprintf(’kd = %f\n’,kd);
fprintf(’kff = %f\n’,kff);%% Reference for ball position
r1 = 20-5;
r2 = 20+5;%% LPF カットオフ周波数
wf = 2*pi*5;%% LPF for ball position sensor

7475
76
77
78
79
80
81
82
83

Fc = wf^2/(s^2+2*0.7*wf*s+wf^2);
Fd = c2d(Fc,ts,’tustin’);[numlpf,denlpf] = tfdata(Fd,’v’);%% Open simulink model
open_system(’bb_pid_hg_sl’);
open_system(’bb_pid_hg_sl/Servo angle’)
open_system(’bb_pid_hg_sl/Ball position and reference’)%% EOF of bb_pid_hg.m

■ 7.3.2

最適状態フィードバックベクトルの設計 (p.190)
プログラム 7-1 [bb lqr test.m]

1234
567
8
9101112
13
1415
16
1718
1920
21
22232425
26
2728
2930
31
32
33
34
35
36
37
38
3940
41
42
43
4445
46
47
48
49
50
5152
5354555657
58596061626364
65%% bb_lqr_test.m%% Initialize & load dataclose allclear allload sim_param%% Set identified parameters
K = 300;
T = 0.28;%% ビーム傾き角 [rad] -> ボール位置 [cm] までの 1/s^2 のゲイン
K_b = (3/5*9.8)*100;%% サーボ 1 度 あたりのアームの傾き% 自作サーボの場合
K_theta_hg = (pi/180)*(2.1/15); % [rad/deg]%% 状態方程式
A = [0 1 00;
0 0 K_b*K_theta_hg 0;
0 0 01;0 0 0-1/T ];B = [ 0;0;
0;
K/T ];
C = [ 1 0 0 0 ];D = 0;
Pbb = ss(A,B,C,D);%% LQ design 1
Q1 = diag([100 20 20 1]);R1 = 10000;
F1 = lqr(Pbb,Q1,R1);%% LQ design 2
Q2 = diag([500 20 20 1]);R2 = 10000;
F2 = lqr(Pbb,Q2,R2);%% Simulationt = 0:ts:5;
x0 = [10,0,0,0]’; % Set initial state% Closed-loop system
Pcl1 = ss(A-B*F1,zeros(4,1),C,D);
Pcl2 = ss(A-B*F2,zeros(4,1),C,D);% Initial state response[y1,tt1,X1] = initial(Pcl1,x0,t);[y2,tt2,X2] = initial(Pcl2,x0,t);%% Plot figure
figure(1)subplot(211)
plot(t,X1(:,1:2),t,X2(:,1:2),’--’)
xlabel(’Time [s]’), ylabel(’Ball position and velocity’)legend(’Design 1 x’,’Design 1 d/dt x’,’Design 2 x’,’Design 2 d/dt x’)
subplot(212)plot(t,X1(:,3:4),t,X2(:,3:4),’--’)
xlabel(’Time [s]’), ylabel(’Servo angle and angular velocity’)legend(’Design 1 \phi’,’Design 1 d/dt \phi’,’Design 2 \phi’,’Design 2 d/dt \phi’)%% EOF of bb_lqr_test.m■ 7.4.4

実験してみよう (p.198–200)
プログラム 7-2 [bb lqr servo hg.m]12345678910111213
1415
16171819202122232425
262728
2930
3132
33
34
35
36
37
38
3940
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

%% bb_lqr_servo_hg.m%% Initialize & load dataclose allclear allload sim_param%% Set identified parameters
K = 300;
T = 0.28;%% Discrete-time plant model
P = K/(T*s^2 + s);
Pd = c2d(P,ts,’zoh’);[numpd,denpd] = tfdata(Pd,’v’);%% ビーム傾き角 [rad] -> ボール位置 [cm] までの 1/s^2 のゲイン
K_b = (3/5*9.8)*100;
Pb = K_b/s^2;
Pbd = c2d(Pb,ts,’zoh’);[numbd,denbd] = tfdata(Pbd,’v’);%% サーボ 1 度 あたりのアームの傾き% 自作サーボの場合
K_theta_hg = (pi/180)*(2.1/15); % [rad/deg]%% 状態方程式
A = [0 1 00;
0 0 K_b*K_theta_hg 0;
0 0 01;0 0 0-1/T ];B = [ 0;0;0;
K/T ];
C = [ 1 0 0 0 ];D = 0;Pbb = ss(A,B,C,D);%% LQI design
Qe = diag([100 20 20 1]);
ru = 10000;
rv = 1000;% Error system
Abar = [A,B;zeros(1,4), 0 ];Bbar = [ zeros(4,1) ;1];
Pbar = ss(Abar,Bbar,[],[]);
Q= [ Qe,zeros(4,1) ;zeros(1,4), ru];
Fbar = lqr(Pbar,Q,rv);
E= [A, B;
C, zeros(1,1) ];Ftmp = Fbar*inv(E);F= Ftmp(:,1:4);
Ki= Ftmp(:,5);Faug = [F, -Ki];%% 閉ループ系の構成とステップ応答% 拡大システムAaug = [ A, zeros(4,1) ;-C, 0];
Baug = [ B ;0 ];Caug = [ C, 0 ];Daug = [ 0 ];Br= [ zeros(4,1) ;1]; % 目標値に対する入力行列

7475
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108

% 閉ループシステム
Paug = ss(Aaug-Baug*Faug,Br,Caug,Daug);% ステップ応答計算
t=0:ts:5;[y,tt,x] = step(Paug,t);%% Plot figure
figure(1)subplot(211)
plot(t,x(:,1:2))
xlabel(’Time [s]’), ylabel(’Ball position and velocity’)legend(’x’,’d/dt x’)
subplot(212)plot(t,x(:,3:4))
xlabel(’Time [s]’), ylabel(’Servo angle and angular velocity’)legend(’\phi’,’d/dt \phi’)%% Reference for ball position
r1 = 20-5;
r2 = 20+5;%% LPF カットオフ周波数
wf = 2*pi*5;%% LPF for ball position sensor
Fc = wf^2/(s^2+2*0.7*wf*s+wf^2);
Fd = c2d(Fc,ts,’tustin’);[numlpf,denlpf] = tfdata(Fd,’v’);%% Open simulink model
open_system(’bb_lqr_servo_hg_sl’);
open_system(’bb_lqr_servo_hg_sl/Ball position and motor angle’)%% EOF of bb_lqr_servo_hg.m