Frequency response data based disturbance observer design
applicable to non-minimum phase systems
Xiaoke Wang 1† , Wataru Ohnishi1 , and Takafumi Koseki1
1

Department of Electrical Engineering, The University of Tokyo, Tokyo, Japan(Tel:+81-3-5841-6676; E-mail:wangxiaoke@koseki.t.u-tokyo.ac.jp;ohnishi@ieee.org;takafumikoseki@ieee.org)
Abstract: Disturbance observer (DOB) has been widely used in industrial field due to it’s simplicity and effectivenessin rejecting disturbance. The performance of disturbance observer is greatly influenced by the bandwidth of low passfilter (Q filter). This paper proposes a novel way for optimizing bandwidth of Q filter by considering experimentallyobtained frequency response data (FRD) of plant. By transforming all the nonliner constraints into convex constraints, theconvex optimisation method can be employed to solve this problem easily. The simulation results verified the feasibilityof proposed method and demonstrated that the optimized Q filter can gurantee the disturbance rejection performance.
Keywords: Disturbance observer, Q filter bandwidth, frequency response data, convex optimisation

1. INTRODUCTION
In industrial applications, unavoidable disturbances
deteriorate the performance of designed control system.
To estimate and compensate the effect of disturbance, disturbance observer has been proposed [1]. As an effective
tool, DOB has been used in various applications, suchas robot manipulators [2], high-speed positioning stage[3] etc. In disturbance observer configuration, Q filter
is necessary to guarantee the causality of system and the
bandwidth of Q filter is desired to be as high as possible
to ensure satisfactory disturbance rejection performancein a wide frequency range. However, the bandwidth of Q
filter is limited by system robustness and noise, and thus
can not be shaped freely [4]. Moreover, in case of a nonminimum phase plant, internal instability due to unstablezeros sets additional limitation for Q filter design [5].
Various guidelines have been proposed on DOB design. By multiplying a filter to shape the plant into nominal plant, Q filter’s bandwidth is increased in [6]. State
space representation of plant is used to design Q filter in[7]. Focusing on closed loop instead of inner DOB loop
to ensure disturbance attenuation performance and robust
stability, Q filter’s design was given in [8]. For nonminimum phase plant, from the prospective of system robustness and sensitivity function limitation, the range of
Q filter’s bandwidth has been given in [5] [9] by using
Bode integral theorem and Poisson integral theorem. In[10], a new filter whose parallel connection with the plant
becomes minimum phase is designed. Then the conventional DOB configuration procedure is employed for the
new system.
The aforementioned methods use the parametricmodel of plant (transfer function or state space representation identified from FRD) which can present the system
dynamics directly. However, since identified parametric model involves model uncertainties, the robustness of
DOB can be deteriorated. Therefore, in this paper, direct usage of frequency response data based DOB systemdesign has been explored.
† Xiaoke Wang is the presenter of this paper.

Previous frequency response data based researchmainly focuses on designing linearly parameterized fixedorder feedback controller by loop shaping method. In[11], authors define margin constraint which is linear withrespect to controllers’ parameters and obtain controller
by linear programming. In [12], non-convex optimization
method is employed to derive controller. Convex optimization has been used to compute robust controllers for
single-input-single-output systems depicted by frequency
response data in [13] which has been applied to specifically design PID controller in [14]-[15], and further extended to multi-input-multi-output case in [16]. Based on
this method, Matlab toolbox has been developed in [17].
In this paper, integrating the frequency response data
based loop shaping method into DOB design has been
researched.
1. DOB system for a non-minimum phase plant isdesigned based on the frequency response data ofplant. The bandwidth of DOB has been optimized
by employing FRD based loop shaping optimization
method.
2. A general way of deriving convex constraints fromoriginal nonlinear constraints which limit the peakvalue of sensitivity function and complementarysensitivity function was given. By introducing conservatism into optimization, this problem can bechanged into convex optimization problem and be
solved easily.
Remaining part of this paper is organized as follows.
Section 2 provides problem formulation, nonlinear constraints will be derived in this section followed by mathematical transformation to convex constraints in Section
3. Based on the constraints obtained, some simulation results are given in Section 4. This paper ends by giving
concluding remarks and future work in Section 5.

2. PROBLEM FORMULATIONIn the disturbance observer control system as shownin Fig. 1, Pr and Pn denote real plant and nominal minimum phase plant, defined by FRD and transfer function, respectively. Q represents the to-be-designed filter.

ˆ r, y are disturbance external input, estimated disturd, d,bance, reference input and output, respectively. The feedback controller Cf b is assumed to be given.

r+ e C
fb
−

d

+

+
−

Q

1

ωp

(4a)

Subject to

0 < ωp < τ1 < ωt ,0.5 < k < 1,
|Wp S| < 1,
|Wm T | < 1.

(4b)
(4c)
(4d)
(4e)

τ,k,ωt

y

−

+

dˆ

Pn−1 QIn this paper, Q is selected as follows in which τ and k
are parameters to be decided and the relative order of Pn
is 2.

τ 2 s2 + 2kτ s + 1

MaximisePr
FRD

Fig. 1 Block diagram of disturbance observer control
systemQ=

In summary, the problem can be formulated into,(1)

then the following equations can be obtained for the blue
dotted part in Fig. 1.
L = Pn−1 Q(1 − Q)−1 Pr (jω)
Pr (jω)Pn−1
N
≜ ,
(2a)
2
2
(s τ + 2kτ s)
D
D
1
≜
,
S=
−1
−1
D+N
1 + (1 − Q) QPn Pr (jω)
(2b)
=

y
Pr (jω)== SPr (jω),
d
1 + (1 − Q)−1 QPn−1 Pr (jω)
(2c)
−1
−1
(1 − Q) QPn Pr (jω)
dˆ=
= 1 − S = T,d
1 + (1 − Q)−1 QPn−1 Pr (jω)(2d)

For nonlinear constraints in terms of optimization variables, in next section, transformation to convex constraints will be given.

3. MATHEMATICAL
TRANSFORMATION OF NONLINEARCONSTRAINTS
In this section, the nonlinear constraints are all transformed into linear function or Linear Matrix Inequality(LMI) form of variables τ, k, ωs , ωt which can simplify
the problem significantly.
3.1. Preliminary
Schur Complement has been used throughout the paper which is introduced as follows where A, B, C, D arescalars or matrices and B ∗ is the conjugate transpose of
B.
[
]
A B
> 0 ⇔ A > 0, A − BD−1 B ∗ > 0. (5)B∗ D
Besides, linear approximation is extensively applied to
the following section. The basic concept is to estimate
the value of a function, f (x), near a point x0 , using thefollowing formula in which f ′ (x0 ) means the slope of thetangent line at x0 .
f (x) ≈ f (x0 ) + f ′ (x0 )(x − x0 )(6)

Additionally, |A| denotes the magnitude of A in which
A is a function.
3.2. Constraint Eq. (4b)1
is a convex funcτition, linear approximation can be directly used. The approximated value is always smaller than the original one,
therefore, if the new inequality holds, the original one
must be satisfied.
2
1
−τi
≥ 2 +
,
(7)
τi
τi−1
τi−1
For the left side of Eq.(4b), sincein which N = Pr (jω)Pn−1 , D = s2 τ 2 + 2kτ s and
L, S, T represent the open loop function, sensitivity
function and complementary sensitivity function, respectively. jω means sequential frequency points.
From above equations, to obtain good disturbance rejection performance, the peak value of sensitivity function should be limited and the 0 dB crossover frequencyof sensitivity function should be as large as possible. By
selecting weighted function for S and T each, the peakvalue of sensitivity and complementary sensitivity function will be no larger than 2 and 1.25 (commonly used inindustrial field). The ωp and ωt will be decided by optimization process.
0.5s + ωp, |Wp S| < 1,ss + ωt, |Wm T | < 1.
Wm =1.25ωt
Wp =

(3a)
(3b)

in which τi means the current value while τi−1 means theprevious value in the optimization loop (τi−1 > τi ).
Therefore, the original constraint can be changed into
following form.
0 < ωp <

2
−τi
+
.
2
τi−1
τi−1

(8)

By using Schur Complement, the right side of Eq. (4b)can be rewritten as the following LMI form.
[
]
1
ωt 1
< ωt ⇔
> 0.
(9)
1 τ
τ

3.3. Constraint Eq.(4d)For the sensitivity constraint, following method is used
to obtain LMI form [18].

In summary, the original nonlinear constraint (4d) canbe transformed into the following LMIs by combining(11), (12), (14), and (16).


0.5s + ωp
D < |D + N (jω)| . (10)s

|Wp S| < 1 ⇔2

2

D + N (jω)
<
s


2

0.5s+ωp
s



⇔
>0
0.5s + ωp ∗
D+N 2
)
(
s
s
[
]
S11
S12
=
> 0.
(S12 )∗ S22
D 2
s

(11)

S11 =
≥

1

1

1

s + ωt1.25ωt(12)⇔

τi−2 still needs to be dealt with. By using the following
technique, the lower bound of τi−2 can be obtained .

(19)

(14)

1
ϕ1 .
2
(ω 2 τ(i−1)
+ 4)

T11
(T12 )∗,

2<

D + N (jω)
sN (jω)s2,(s + ωt ) > 0,2 
D + N (jω) 
1.25

1.25
]
T12
> 0.
T22

s

As before, T11 and T22 needs transformation. Since
T22 = S22 , this part is omitted due to the repeatence. For
T11 , the linear approximation is used.

(15)

As for S22 , the linear approximation was employed.
Let N (jω)= xk + yk j,s2

S22(i) = (sτi2 + 2ki τi ) + (xk + yk j)′(τi − τi−1 )= S22(i−1) + S22(τ,i−1)(16)′(ki − ki−1 ) = Φ,+ S22(k,i−1)

in which22+ 2ki−1 τi−1 ) + (xk + yk j) ,S22(i−1) = (sτi−1
(17a)
2
′
= 4ωτi−1 (τi−1ω + yk )S22(τ,i−1)+ 4ki−1 (2ki−1 τi−1 + xk ),′S22(k,i−1) = 4τi−1 (2ki−1 τi−1 + xk ).

=

2

(20)

In conclusion. S11 part can be transformed into≥

D + N (jω)
N (jω)

2 N (jω)s2|1.25ωt |[−2
By introducing a new variable ϕ1 and making 2τi−1−−4 2
τi−1 τi > ϕ1 > 0, LMI form constraint can be obtained.

]τi> 0, ϕ1 > 0.
1

<

2

|ωt |


⇔  N (jω)(s + ωt ))∗
( s

−2 −2
−4
−2
−4 2
⇔ τi−4 ≥ 2τi−1
τi − τi−1⇔ τi−2 ≥ 2τi−1− τi−1
τi .
(13)

[ 2
4
2τi−1 − ϕ1 τi−1
τi2

|s + ωt |−2−2(τi−2 − τi−1)(τi−2 − τi−1)≥01

s + ωtD + N (jω)<.
1.25ωtN (jω)

After squaring both sides and using Schur Complement,
the following results can be developed.

.
= 2 2
2
ω 2 τi2 τ(i−1)(ω τ(i−1) + 4) τi2+ 4τi2

D 2
s

(18b)

1

1

2 = ω 2 τ 4 + 4k 2 τ 2|(sτ 2 + 2kτ )|

S11 =

(18a)

3.4. Constraint Eq. (4e)For complementary sensitivity constraint, by following similar process as dealing with previous one, the desired LMI form constraint can be obtained.
|Wm T | < 1 ⇔S11 and S22 still need to be transformed into linear
function of variables. For S11 ,0.5s+ωpsϕ12 (ω 2 τ(i−1)+ 4) > 0,
0.5s + ωp ∗
(
)
Φ
s
[ 2
]
4
2τi−1 − ϕ1 τi−1
τi> 0, ϕ1 > 0.
τi1

Square on (10) both sides and turn this inequality into
matrix inequality form by using Schur Complement.
0.5s + ωp
D
s
s

1



1

(17b)
(17c)

2
2
ωt(i)
≥ 2ωt(i−1) ωt(i) − ωt(i−1),(21)

in which ωt(i) is the current value while ωt(i−1) represents for the previous one in optimization loop.
Combining (16), (20) and (21), the original nonlinear
constraint (4e) can be changed into
2ωt(i−1) ωt(i) − ω 2t(i−1)
 N (jω)
∗
(
(s
+
ω
t(i) ))s
1.25

N (jω)s(s + ωt(i) )
1.25




 > 0.
Φ
(22)

3.5. Constraints summary
After finishing all the process mentioned above, theoriginal constraints have been formulated into a new form

as follows.

1

ωp

0.5

2
i
Subject to ωp > 0, ϕ1 > 0, ωp < τ−τ+ τi−1, (23b)
2
i−1
[
]
2
4
0.5 < k < 1, 2τi−1 − ϕ1 τi−1 τi > 0,
(23c)
τi
1

1
0.5s + ωp 
ϕ
[
]
1
2
 (ω 2 τ(i−1)

ωt 1
+ 4)
s
 > 0,(23d)1 τ > 0,  0.5s + ωp
∗
(
)
Φ
s


N (jω)
(s
+
ω
)
t(i)
s
2ωt(i−1) ωt(i) − ω 2
t(i−1)

1.25
 N (jω) > 0.(23e)∗((s + ωt(i) ))s
Φ
1.25
By finishing the reformulation of the problem, theoptimization problem has been changed in to a convex
optimization problem and can be solved by commercial
solvers, for instance yalmip plus Mosek[19][20] in Matlab.

0

-2.5
-1.5

Pr (jω)(s + 2.101).
2
84.977(sτ + 2kτ ) + Pr (jω)(s + 2.101)

Magnitude [dB]

(24d)

102

101

102

Phase [deg]

270
180
90
0
-90
10-1

100

-0.5

0

0.5

1

Magnitude [dB]Fig. 3 Nyquist plot of before and after optimization
20
0
-20
-40
10-1

100

101

102

101

102

90
0
-90
-180
100

Frequency [Hz]Fig. 2 Bode plot of real plant(frd) and nominal plant

Magnitude [dB]Fig. 4 Bode plot of open loop function L, Wp S and
Wm T before optimization
20
0
-20
-40
10-1

100

101

102

101

102

90

Phase [deg]84.977
1,Q = 2 2, (24a)s(s + 2.101)
τ s + 2kτ s + 1(s + 2.101)
L=Pr (jω),(24b)
84.977(sτ 2 + 2kτ )84.977(sτ 2 + 2kτ )S=,
84.977(sτ 2 + 2kτ ) + Pr (jω)(s + 2.101)
(24c)

101

-1

Frequency [Hz]

Pn =

100

Before optimisation
After optimisation

Real axis

-270
10-1

4.1. Simulation plant
In this paper, the Bode plot of nominal plantPn (2nd order transfer function) and actual real plant
Pr (jω)(frequency response data) is shown in Fig. 2. The
nominal plant catches the dynamics of real plant in lessthan 10Hz frequency range and doesn’t contain unstablezeros.

40
20
0
-20
-40
-60
-80
-100
10-1

-1

-2

4. SIMULATION RESULT AND
ANALYSIS

T =

-0.5

-1.5

Phase [deg]τ,k,ωt ,ϕ1

Imaginary axis

(23a)

Maximise

0
-90
-180
-270
10-1

100

Frequency [Hz]Fig. 5 Bode plot of open loop function L, Wp S and
Wm T after optimization
4.2. Case study result
In this section, by building the constraints following
the proposed method, the ωp can be maximized. The
simulation initial condition is given as τinit = 0.1 [s],kinit = 0.9, ωp(init) = π [rad/s], ωt(init) = 80π [rad/s].
After using convex optimization, τopt = 0.0419 [s],kopt = 1, ωp(opt) = 9.9 [rad/s]. In Fig. 3, dashed blackline represents for unit circle and dotted black line is a
circle whose center locates at (-0.5,0) and the radius is0.5. Since the closest distance from Nyquist plot to critical point (−1, 0) is modulus margin (the inverse of the
maximum of the sensitivity function), Nyquist plot has nointersection with dashed circle implies that the sensitivity function constraint holds successfully. And proposed
optimization method makes the Nyquist plot become tangent to the dashed circle. The predefined peak value con-

0.5

0

0.4

-20
-40
25

30

Output response(mm)Magnitude [dB]20

35

Phase [deg]

90
0
-90
-180

No DOB
InitialOptimized

0.3
0.2
0.1
0
-0.1

-270
25

30

35

-0.2

Frequency [Hz]

0

0.5

1

1.5

2

2.5

3

3.5

4

Time(s)Fig. 6 Locally enlarged figure of Fig. 5

Fig. 9 Output response when initial Q filter and optimized Q filter are employedstraints for S and T have been satisfied can also be told
more straightforwardly from Fig. 5 and Fig. 6 in which|Wp S| and |Wm T | are always under 0 dB.
4.3. Disturbance rejection performance
Simulations also have been done to test the disturbancerejection performance after using designed Q filter in Fig.
7. During this simulation, a well-identified transfer function (tf) of Pr (8th order) (as shown in Fig. 8) was used
and feedback controller Cf b was chosen as follows.

r+ e C
fb
−

+

+

d

−

Q

y
Pr
tf
−

+

dˆ

Pn−1 Q

Magnitude [dB]Fig. 7 Block diagram of disturbance observer control
system
40
20
0
-20
-40
-60
-80
-100
10-1

5.1. Conclusion
This paper has proposed a novel optimization method
for maximizing bandwidth of disturbance observer configuration by just using frequency response data. What’smore, all the nonlinear constraints are transformed into
convex form which can be solved by convex optimization method. The simulation results verify the feasibilityof proposed method.
5.2. Analysis and Future work
In this paper, two parameters of Q filter has been
tuned. Although the real plant is a non-minimum phaseone, the nominal plant in this paper is a second order minimum phase one. In the future, proposed method will
be extended to the design of DOB system when nonminimum phase nominal plant is involved.

6. ACKNOWLEDGMENT

100

101

102

100

101

102

The work is supported by JSPS KAKENHI Grant
Number 18H05902 and the Telecommunications Advancement Foundation. The student is supported by
China Scholarship Council (CSC). The authors thankProfessor Hiroshi Fujimoto (The University of Tokyo)
cordially for providing the frequency data.

270

Phase [deg]5. CONCLUDING REMARKS

180
90
0
-90
10-1

REFERENCES[1]Frequency [Hz]Fig. 8 Bode plot of Pr (8th order transfer function) andPr (frequency response data)[2]
Cf b = 1.36 +3.64
0.261s+.
s0.0249s + 1(25)

Making reference input as zero and injecting unit step disturbance to the system, the output response are shown in
Fig. 9. The proposed Q filter design outperforms initial
Q filter and disturbance rejection performance has been
improved.

[3]

K. Ohishi, K. Ohnishi, and K. Miyachi, “Torquespeed regulation of DC motor based on load torqueestimation method ”, International Power Electronics Conference, pp. 1209–1218, 1983.
E. Schrijver and J. van Dijk, “Disturbance observers
for rigid mechanical systems: equivalence, stability, and design”, Journal of Dynamic Systems, Measurement, and Control, vol. 124, no. 4, pp. 539-548,
2002.
S. Endo, H. Kobayashi, C. J. Kempf, S. Kobayashi,
M. Tomizuka, and Y. Hori, “Robust digital tracking controller design for high-speed positioning systems”, Control Enigineering Practice, vol. 4, no. 4,
pp. 527–536, 1996.

[4]K. Ohnishi, M. Shibata, and T. Murakami, “Motion
control for advanced mechatronics”, IEEE/ASME
Transactions on Mechatronics, vol. 1, no. 1, pp. 56–
67, 1996.
[5] E. Sariyildiz and K. Ohnishi, “A new solution for
the robust control problem of non-minimum phasesystems using disturbance observer”, IEEE International Conference on Mechatronics, pp. 46–51,
2013.
[6] X. Chen and M. Tomizuka, “Optimal plant shaping
for high bandwidth disturbance rejection in discrete
disturbance observers”, Proceedings of the American Control Conference, pp. 2641–2646, 2010.
[7] H. Shim and Y. J. Joo, “State space analysis of disturbance observer and a robust stability condition”,
Proceedings of the IEEE Conference on Decision
and Control, pp. 2193–2198, 2007.
[8] C. C. Wang and M. Tomizuka, “Design of robustly
stable disturbance observers based on closed loop
consideration using H ∞ optimization and its applications to motion control systems”, Proceedings of
the American Control Conference, pp. 3764–3769,
2004.
[9] E. Sariyildiz and K. Ohnishi, “A Guide to design
disturbance observer”, Journal of Dynamic Systems,
Measurement, and Control, vol. 136, no. 2, 2014.
[10] H. Shim, N. H. Jo and Y. Son, “A new disturbanceobserver for non-minimum phase linear systems”,
Proceedings of the American Control Conference,
pp. 3385-3389, 2008.
[11] A. Karimi, M. Kunze, and R. Longchamp, “Robust
controller design by linear programming with application to a double-axis positioning system”, Control
Engineering Practice, vol. 15, no. 2, pp. 197–208,
2007.
[12] A. J. Den Hamer, S. Weiland, and M. Steinbuch,“Model-free norm-based fixed structure controller
synthesis”, Proceedings of the IEEE Conference on
Decision and Control, pp. 4030–4035, 2009.
[13] A. Karimi and G. Galdos, “Fixed-order H ∞ controller design for nonparametric models by convex
optimization”, Automatica, vol. 46, no. 8, pp. 1388–
1394, 2010.
[14] M. Hast, K. J. Åström, B. Bernhardsson, and
S. Boyd, “PID design by convex-concave optimization ”, Proceedings of the European Control Conference, pp. 4460–4465, 2013.
[15] K. Nakamura, K. Yubai, D. Yashiro, and S. Komada, “Controller design method achieving maximization of control bandwidth by using Nyquist
diagram ”, Proceedings of the International Automatic Control Conference, pp. 35–40, 2016.
[16] G. Galdos, A. Karimi, and R. Longchamp, “H ∞controller design for spectral MIMO models by
convex optimization” , Journal of Process Control,
vol. 20, no. 10, pp. 1175-1182, 2010.
[17] A. Karimi, “Frequency-domain robust control toolbox ” , Proceedings of the IEEE Conference on De-

cision and Control, pp. 3744-3749, 2013.
[18] K. Nakamura, K. Yubai, D. Yashiro and S. Komada,“Fully parameterized controller design method forhigh control bandwidth using frequency response
data sets ”, Proceedings of the International Workshop on Sensing Actuation, and Motion Control,
2017.
[19] J. Löfberg, “YALMIP : a toolbox for modelingand optimization in MATLAB”, IEEE International
Conference on Robotics and Automation, pp. 284289, 2004[20] MOSEK ApS, “MOSEK optimization toolbox for
MATLAB manual”, 2019.