• 대한전기학회
Mobile QR Code QR CODE : The Transactions of the Korean Institute of Electrical Engineers
  • COPE
  • kcse
  • 한국과학기술단체총연합회
  • 한국학술지인용색인
  • Scopus
  • crossref
  • orcid
The Transactions of the Korean Institute of Electrical Engineers

The Transactions of the Korean Institute of Electrical Engineers

ISO Journal TitleTrans. Korean. Inst. Elect. Eng.
SCImago Journal & Country Rank

Journal Search

  1. Home
  2. Archive
  3. 2018-08 (Vol.67 No.08)
  4. 10.5370/KIEE.2018.67.8.1062
XML PDF INFO REF
[

Research article

]

The Transactions of the Korean Institute of Electrical Engineers

KIEE Vol. 67, No. 8, p.1062-1070

ISSN (print) :

1975-8359

ISSN (online) :

2287-4364

Received : 21 July 2018Revised : Accepted : 30 July 2018

DOI :

http://doi.org/10.5370/KIEE.2018.67.8.1062

펄스응답에 의한 저차 이산시간 모델의 식별

Identification of Discrete-Time Low-Order Model from Pulse Response

황지호 (Jiho Hwang) 1 차승표 (Seungpyo Cha) 1 김영철 (Young-Chol Kim) 2†
  1. ( Dept. of Electronics Eng., Graduate School, Chungbuk National University, Korea )
  2. ( Dept. of Electronics Eng., Chungbuk National University, Korea )

† Corresponding author, E-mail: yckim@cbu.ac.kr

Copyright © The Korean Institute of Electrical Engineers(KIEE)

License :

This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.( www.kiee.or.kr ).

Abstract

This paper presents a simple identification method for discrete-time low-order model of unknown delay process from pulse response. The key idea is to find the parameters of the model such that the first N moments of the unknown process and the model are equal. We first show that the k-th moment of a process can be determined by the moments of the input and output. The parameters and delay are estimated separately. It is shown that for a given delay, the parameters of the low-order model can be determined by solving linear equations in a matrix form. Delay of the model is estimated such that the integral of the absolute errors (IAE) of the candidate models with possible delays minimizes. The illustrative example shows that the proposed method can directly identify low-order models without order reduction process from a single pulse response.

Key words

Method of moment matching, Parameters identification, Discrete-time low-order models

1. 서론

정교한 성능을 낳는 뛰어난 현대적 제어기법이 개발되어있음에도 불구하고 산업용 프로세스 제어기는 95% 이상이 PID와 1차 보상기 형태의 저차 제어기를 사용하고 있다(1). 다시 말해, 이러한 프로세스는 대부분 2차 이하의 저차 모델로 표현될 수 있음을 반증하는 것이라 할 수 있다. 또한 시스템이 원하는 제어 성능을 만족하도록 제어기를 설계하고 이를 분석하기 위해서는 플랜트의 수학적 모델이 필요한데 특히, 디지털 제어기를 설계할 경우 플랜트에 대한 이산시간모델을 구해야 한다. 이러한 관점에서 프로세스를 2차이하의 이산시간 전달함수 모델로 식별하는 것은 매우 흔히 나타나는 일이다. 여기서는 플랜트에 간단한 펄스입력을 인가하고 그 응답으로부터 이산시간 저차 전달함수 모델을 직접 식별하는 문제에 주목하고자 한다.

대표적인 입력-출력 이산 시간 모델링 방법은 소위 최적화 방법(optimization methods)으로서 최소자승추정(least square estimation: LSE), 기구변수방법(instrumental variable method: IV), 예측오차방법(prediction error method: PEM) 등(2,3)이 있다. 가장 정밀한 추정 모델을 제공하지만 알고리즘의 구현이 복잡하고 시험입력은 지속여기조건(persistent excitation)을 만족해야 한다. 더욱이 이러한 최적화 기법들은 샘플링 순간에서 입-출력 데이터를 모델 수식에 맞추는(fitting) 방식이기 때문에 식별 모델의 차수가 작게 선택되거나 출력 잡음 정보(확률분포, 백색성)가 충분치 않으면 무시할 수 없는 모델링 오차를 유발할 수 있다. 이러한 이유로 식별된 모델이 불필요하게 고차일 수도 있고, 또 저차 제어기 설계 시에는 플랜트 모델의 차수 줄임(model order reduction) 과정(4-6)이 요구된다. 차수 줄임 기법을 이용하여 저차 모델로 표현할 때, 일반적으로 저차 모델의 지연시간을 고차 모델의 지연시간과 동일하게 놓는데 이 지연시간은 저차 모델을 위한 최적의 지연시간이 아닐 수도 있음을 유의해야 한다.

더욱이 연속시간 프로세스의 지연시간은 일반적으로 샘플링 시간의 정수배가 되지 않으며 분수지연시간 (fractional time- delay)은 이산시간모델의 계수에 나타나기 때문에 모델 파라미터와 샘플링 시간, 지연시간은 서로 종속되어 있다. 다음 연속시간 1차 지연시간 프로세스의 이산시간 모델을 예로 들어보기로 한다.

연속시간 FOTD모델: G ( s ) = K e - L s τ s + 1

샘플링 시간 Ts, 지연시간은 L = d T s + λ T s ,   ( 0 ≤ λ ≤ 1 ) 이라 하자 (여기서 d는 정수). 위 모델의 이산시간모델을 ZOH (zero-order hold)를 포함한 펄스전달함수로 표현하면

이산시간 FOTD모델: H ( z - 1 ) = b 1 z - 1 + b 2 z - 2 1 + a 1 z - 1 z - d

여기서,

b 1 = K ( 1 - e ( λ - 1 ) T s τ ) ,     b 2 = K ( e ( λ - 1 ) T s τ - e - T s τ ) ,     a 1 = - e - T s τ

위 식으로부터 λ = 1 ,   ( L = ( d + 1 ) T s ) 이면 b1=0이고, H(z-1)의 지연시간은 만큼 Ts증가한다. 또, λ=0이면 b2=0가됨을 알 수 있다. 분수 지연시간이 존재하면 분자의 계수가 변하고 이산시간 모델의 영점이 변하게 된다. 이 예는 고차인 실제 프로세스를 직접 저차모델로 식별하려 할 때 정수 지연시간 d와 파라미터 {b1 b2 a1}는 상호 결합되어 있기 때문에 이 둘을 동시에 추정하는 하는 것은 매우 어려운 문제임을 드러낸다.

최근 Kim 등(7)은 모멘트정합 기법(method of moment matching)을 이용하여 연속 시간 저차 모델식별 알고리즘을 제시하였다. 모멘트 정합은 주로 고차 모델을 저차 모델로 근사화하기 위한 유력한 방법으로 사용되어 왔다(5,6). 이는 저차 모델의 0~k차 모멘트를 고차 모델의 0~k차 모멘트와 서로 정합 되도록 저차 모델의 파라미터를 결정하는 방식이다. Kim 등은 이 개념을 확장하여 입출력 응답의 모멘트로부터 직접 연속시간 저차모델의 파라미터를 식별하는 공식을 유도하였다. 이 방법은 샘플 시각에서의 입출력 데이터를 수식에 맞추는 최적화 방법에 비해 시험 입력이 구형 펄스처럼 단순한 입력이어도 매우 정확한 저차 모델을 직접 얻을 수 있음을 보였다. 이 방법의 가장 큰 장점은 모멘트 정합에 의한 차수 줄임을 내재적으로 만족시키며 모델 파라미터를 식별한다는 것이고 또한 시험입력의 지속여기조건을 요구하지 않는다는 것이다. 이종언 등(8)은 자율 주행 자동차(투싼ix, 2000cc, 가솔린, 전륜구동, 2010, 현대자동차)의 가감속 모델을 구하는데 이 방법을 적용하여 실제 문제에 쉽게 적용할 수 있음을 보였다.

본 논문의 주요 목적은 이러한 모멘트 정합 방식을 이산 시간 저차 모델을 식별하는 문제에 적용하여 간단한 펄스응답으로부터 직접 이산시간 저차모델을 식별하는 알고리즘을 제시하는데 있다. 또한 이결과로부터 가장 많이 사용되는 1차지연시간(first- order plus time-delay: FOTD), 2차지연시간(second-order plus time-delay: SOTD), 2차무지연시간(second-order with a free of delay: SODF) 모델 파라미터의 식별 공식을 제시한다. 연속 시간 저차 모델 식별을 위한 모멘트 정합 기법의 장점을 그대로 유지할 수 있기 때문에 매우 유용할 것으로 기대된다.

논문의 구성은, 2장에서 이산시간 신호의 모멘트를 정의하고 전달함수의 도함수와 시스템 모멘트 사이의 관계를 유도한다. 또한 프로세스의 모멘트를 입출력 모멘트로부터 결정할 수 있음을 유도한다. 3장에서는 이산시간 저차모델의 파라미터를 식별하는 공식을 유도한다. 4장에서 적용 수치 예를 보이고 5장에서 결론을 맺는다.

2. 모멘트의 정의와 예비결과

2.1 이산시간 시스템의 모멘트

이산시간 시퀀스 $f[i]$의 k차 모멘트를 다음과 같이 정의한다.

(1)

$M _ { k } := \displaystyle\sum _ { i=0 } ^ { \infty } i ^ { k } f[i], k=1,2,3, \cdots$

단일 입출력 선형 이산시간 시스템 $H(z ^ { -1 } )$를 고려한다. 단, 이 시스템은 최소위상계이고 안정한 시스템이라고 가정한다. 정의 (1)을 이용하면 시스템 $H(z ^ { -1 } )$의 입력 시퀀스 $u(i)$와 출력 시퀀스 $y(i)$에 대한 k차 모멘트는 다음과 같이 나타낼 수 있다.

(2)

$M _ { k } ^ { u } = \displaystyle\sum _ { i=0 } ^ { \infty } i ^ { k } u(i)$, $M _ { k } ^ { y } = \displaystyle\sum _ { i=0 } ^ { \infty } i ^ { k } y(i)$.

전달함수 $H(z ^ { -1 } )$의 임펄스 응답을 $h[i]$라 하면 다음 식으로 표현할 수 있다.

(3)

$H(z ^ { -1 } )= \frac{ b _ { 0 } +b _ { 1 } z ^ { -1 } +b _ { 2 } z ^ { -2 } + \cdots +b _ { m } z ^ { -m } } { 1+a _ { 1 } z ^ { -1 } +a _ { 2 } z ^ { -2 } + \cdots +a _ { n } z ^ { -n } } (m \leq n) $

$= \displaystyle\sum _ { i=0 } ^ { \infty } h _ { i } z^{-i}$.

정의 (1)을 (3)에 적용하면, 프로세스 $H(z ^ { -1 } )$의 k차 모멘트는

(4)

$M _ { k } ^ { h } = \sum _ { i=0 } ^ { \infty } i ^ { k } h_i$.

$z=1$ 부근에서 $H(z ^ { -1 } )$의 $z^ { -1 }$에 대한 k차 도함수를 다음과 같이 정의한다.

(5)

$H ^ { (k) } (z ^ { -1 } ):= \left . \left ( z ^ { -1 } \frac{ d } { dz ^ { -1 } } \right ) ^ { k } H(z ^ { -1 } ) \right | _ { z=1 } ,~ k=0,1,2, \cdots$

이제 이산시간 전달함수의 도함수와 시스템 모멘트와의 관계를 유도하기로 한다. 식 (5)를 $k=0,1,2,3,4$에 대해 전개하고 (4)를 이용하면 다음 관계를 얻는다.

(6)

$H(1)= \displaystyle\sum _ { i=0 } ^ { \infty } h_i =M_0^h$

$H ^ { ' } (1)=z ^ { -1 } \frac{ dH(z ^ { -1 } ) } { dz ^ { -1 } } \biggr | _ { z=1 } = \displaystyle\sum _ { i=0 } ^ { \infty } ih _ { i } z ^ { -i } \biggr | _ { z=1 } =M_1^h $

$H '' (1)= \left . \left ( z ^ { -1 } \frac{ d } { dz ^ { -1 } } \right ) ^ { 2 } H(z ^ { -1 } ) \right | _ { z=1 } = \displaystyle\sum _ { i=0 } ^ { \infty } i(i-1)h _ { i } z ^ { -i } \biggr | _ { z=1 } = M_2^h - M_1^h$

$H ''' (1)= \left . \left ( z ^ { -1 } \frac{ d } { dz ^ { -1 } } \right ) ^ { 3 } H(z ^ { -1 } ) \right | _ { z=1 } = \displaystyle \sum _ { i=0 } ^ { \infty } i(i-1)(i-2)h _ { i } z ^ { -i } \biggr| _ { z=1 }$ $= M_3^h -3M_2^h +2 M_1^h$

$H^ { (4) } (1)= \left . \left ( z ^ { -1 } \frac{ d } { dz ^ { -1 } } \right ) ^ { 4 } H(z ^ { -1 } ) \right | _ { z=1 } = \displaystyle \sum _ { i=0 } ^ { \infty } i(i-1)(i-2)(i-3)h _ { i } z ^ { -i } \biggr | _ { z=1 }$ $= M_4^h -6M_3^h +11M_2^h -6 M_1^h$

식 (6)으로부터 역으로 모멘트를 도함수의 항으로 나타내면

(7)

$M _ { 0 } ^ { h } =H(1)$

$M _ { 1 } ^ { h } =H ^ { ' } (1)$

$M _ { 2 } ^ { h } =H ^ { '' } (1)+H ^ { ' } (1)$

$M _ { 3 } ^ { h } =H ^ { (3) } (1)+3H ^ { '' } (1)+H ^ { ' } (1)$

$M _ { 4 } ^ { h } =H ^ { (4) } (1)+ 6H ^ { (3) } (1)+7H ^ { '' } (1)+H ^ { ' } (1)$.

4차 이상의 도함수에 대해서도 (6), (7)과 같은 방법으로 유도할 수 있으며 본 논문에서 고려하는 SOTD 모델 식별에는 4차 모멘트까지 요구되기 때문에 4차까지 전개하였다. 여기서, 0차모멘트 $M_0^h$는 프로세스 $H(z ^ { -1 } )$의 DC이득임을 알 수 있다.

다음 정리는 프로세스의 k차 모멘트를 입력과 출력의 응답의 k차 모멘트의 관계를 나타낸다.

정리 1: 선형시스템 $H(z ^ { -1 } )$에 유한크기 입력 를 인가하여 출력 를 얻었다고 하자. 시스템 모멘트와 입출력 모멘트 사이에 다음 관계가 성립한다.

(8)

$M _ { k } ^ { y } = \displaystyle \sum _ { i=0 } ^ { k } \begin{pmatrix}k\\i\end{pmatrix} M _ { k-i } ^ { u } M _ { i } ^ { h }$

여기서 $\begin{pmatrix}k\\i\end{pmatrix} = \frac{ k! } { i!(k-i)! }$ (이항 계수)

증명: 선형시스템이므로 입력 에 대한 응답은 컨볼루션합 공식에 의해

(9)

$y(i)= \displaystyle \sum _ { j=0 } ^ { \infty } u(j)h(i-j)$.

식 (9)를 (2)에 대입하면 출력의 모멘트는

(10)

$M _ { k } ^ { y } = \displaystyle \sum _ { i=0 } ^ { \infty } (iT) ^ { k } y(iT)= \displaystyle \sum _ { i=0 } ^ { \infty } (i) ^ { k } \displaystyle \sum _ { j=0 } ^ { \infty } u(j)h(i-j)$

$= \displaystyle\sum _ { j=0 } ^ { \infty } u(j) \displaystyle\sum _ { i=0 } ^ { \infty } (i) ^ { k } h(i-j)$.

여기서 $l=i-j$ 로 치환하고 (10)을 다시 정리하면

$M _ { k } ^ { y } = \displaystyle \sum _ { j=0 } ^ { \infty } u(j) \sum _ { l=0 } ^ { \infty } (l+j) ^ { k } h(l)= \displaystyle \sum _ { j=0 } ^ { \infty } u(j) \displaystyle \sum _ { l=0 } ^ { \infty } \sum _ { i=0 } ^ { \infty } \begin{pmatrix}k\\i\end{pmatrix} (l ^ { k-i } j ^ { i } ) h(l)$

$= \begin{pmatrix}k\\0\end{pmatrix} \displaystyle \sum _ { j=0 } ^ { \infty } u(j) \displaystyle \sum _ { l=0 } ^ { \infty } (l) ^ { k } h(l)+ \begin{pmatrix}k\\1\end{pmatrix} \displaystyle \sum _ { j=0 } ^ { \infty } j u(j) \displaystyle \sum _ { l=0 } ^ { \infty } (l) ^ { k-1 } h(l)+ \cdots $

$+ \begin{pmatrix}k\\k-1\end{pmatrix} \displaystyle \sum _ { j=0 } ^ { \infty } (j) ^ { k-1 } u(j) \displaystyle \sum _ { l=0 } ^ { \infty } l h(l)+ \begin{pmatrix}k\\k\end{pmatrix} \displaystyle \sum _ { j=0 } ^ { \infty } (j) ^ { k } u(j) \displaystyle \sum _ { l=0 } ^ { \infty } h(l)$

$= \displaystyle \sum _ { i=0 } ^ { \infty } \begin{pmatrix}k\\i\end{pmatrix} M _ { i } ^ { u } M _ { k-i } ^ { h }$ ▼

다음 따름정리 1은 프로세스 $H(z ^ { -1 } )$의 k차 모멘트를 입력과 출력의 모멘트로부터 결정할 수 있음을 나타낸다. 증명은 생략한다.

따름정리 1: 이산시간 전달함수 $H(z ^ { -1 } )$에 대해 입력과 출력의 k차 모멘트를 각각 $M _ { k } ^ { u }$, $M _ { k } ^ { y }$ 라 하면 다음 관계가 성립한다.

(11)

$M _ { k } ^ { h } = \frac{ 1 } { M _ { 0 } ^ { u } } \left [ M _ { k } ^ { y } - \displaystyle \sum _ { i=0 } ^ { k } \begin{pmatrix}k\\i\end{pmatrix} M _ { k-1 } ^ { h } M _ { i } ^ { u } \right ] , ~ k=1,2,3, \cdots$

여기서, $M _ { 0 } ^ { h } = \frac{ M _ { 0 } ^ { y } } { M _ { 0 } ^ { u } } .$ ▼

2.2 저차 이산시간 모델 식별 문제

미지 선형 이산시간 프로세스 $H(z ^ { -1 } )$ 에 대해 저차 이산시간 전달함수 모델의 파라미터를 식별하는 문제를 기술한다. 만약 연속시간 프로세스를 고려한다면 < 그림 1>처럼 디지털 입력신호가 ZOH를 통해 인가되고 연속시간 출력을 샘플링 하여 이산시간 입출력 간의 전달함수를 $H(z ^ { -1 } )$라고 간주한다. 앞에서 언급한 바와 같이 프로세스는 안정하고 최소위상계라고 가정한다. 또한 2차이하의 모델로 표현할 수 있는 대상으로 한정한다. 이러한 조건에서 본 논문에서는 간단한 펄스응답으로부터 임의의 이산시간 저차모델 ${ \widehat H } _L (z ^ { -1 } )$의 파라미터와 지연시간을 추정하는 알고리즘을 제시하고 이 결과로부터 다음 FOTD, SOTD 모델을 식별하는 공식을 제시하고자 한다.

(12)

$FOTD : { \widehat H } _ { FD } (z ^ { -1 } )= \frac{ \widehat { b _ { 1 } } z ^ { -1 } + \widehat { b _ { 2 } } z ^ { -2 } } { 1+ \widehat { a _ { 1 } } z ^ { -1 } } z ^ { - \widehat { d } }$

(13)

$SOTD : { \widehat H } _ { SD } (z ^ { -1 } )= \frac{ \widehat { b } _ { 1 } z ^ { -1 } + \widehat { b } _ { 2 } z ^ { -2 } + \widehat { b } _ { 3 } z ^ { -3 } } { 1+ \widehat { a } _ { 1 } z ^ { -1 } + \widehat { a } _ { 2 } z ^ { -2 } } z ^ { - \widehat { d } }$ ${ \widehat H } _L (z ^ { -1 } )= \frac{ N_L (z ^ { -1 } ) } { D_L (z ^ { -1 } ) } z ^ { - \widehat { d } }$

(14)

$= \frac{ \widehat { b } _ { 1 } z ^ { -1 } + \widehat { b } _ { 2 } z ^ { -2 } + \cdots + \widehat { b } _ { q } z ^ { -q } } { 1+ \widehat { a } _ { 1 } z ^ { -1 } + \widehat { a } _ { 2 } z ^ { -2 } + \cdots + \widehat { a } _ { l } z ^ { -l } } z ^ { - \widehat { d } }$

여기서, $\widehat { d }$는 추정할 정수 지연시간이다.

그림 1 이산시간 프로세스 모델

Fig. 1 Discrete-time process model

../../Resources/kiee/KIEE.2018.67.8.1062/fig1.png

2.3 접근 방법

모멘트 정합에 의한 저차모델 식별의 기본개념은 미지 플랜트 $H(z ^ { -1 } )$의 처음 k차 까지 모멘트와 식별하고자 하는 저차 모델 ${ \widehat H } _L (z ^ { -1 } )$의 처음 k차까지 모멘트가 서로 동일한 값을 갖도록 ${ \widehat H } _L (z ^ { -1 } )$를 결정하는 것이다. 이를 구현하기 위해서는 먼저 입출력 데이터로부터 미지 프로세스의 k차 까지 모멘트를 구해야 하는데 이미 앞 절에서 예비결과로 제시한 따름정리 1을 이용하면 된다. 이 때 인가하는 시험입력은 구형펄스, 반주기 정현파처럼 간단한 입력파형이어도 충분하다.

이제 남은 문제는 모멘트 정합을 이루도록 저차모델 ${ \widehat H } _L (z ^ { -1 } )$의 파라미터를 구하는 관계식을 유도하는 일이다. 서론에서 토의한 바와 같이 샘플링 시간과 분수 지연시간의 영향으로 이산시간 모델의 정수 지연시간을 파라미터와 동시에 추정하는 것은 매우 어려운 문제이다. 더욱이 (14) 모델의 $\widehat { d }$는 정수이어야 한다. 여기서는 지연시간과 파라미터를 분리하여 추정하고자 한다. 시간응답으로부터 정수 지연시간의 가능한 후보 값을 선택하는 것은 어렵지 않다. 이 집합을 $D:= \{ d | d_ { j-2 } , d_ { j-1 } , d_ { j } , d_ { j+1 } , d_ { j+2 } , \}$라고 하자. 이 집합의 각 지연시간을 참값으로 간주하여 대응하는 저차모델의 파라미터를 결정한다. 지연시간 후보 집합의 숫자만큼 모델 집합이 얻어지게 된다. 이 모델 중에서 실제 펄스응답과의 누적절대오차(integral of the absolute errors: IAE)를 비교하여 최소 IAE를 주는 지연시간과 파라미터를 최종 식별모델로 결정하게 된다.

3. 주요 결과

본 절에서는 미지 선형 이산시간 시스템의 펄스응답으로부터 2.2절에서 기술한 저차 모델의 파라미터를 결정하는 공식을 유도한다.

따름정리 1에 따라 미지 프로세스 $H(z ^ { -1 } )$의 모멘트 $M _ { k } ^ { h }$는 다음과 같이 입력의 k차 모멘트 $M _ { k } ^ { u }$와 출력의 k차 모멘트 $M _ { k } ^ { y }$의 항으로 표현된다.

(15)

$M _ { 0 } ^ { h } = \frac{ M _ { 0 } ^ { y } } { M _ { 0 } ^ { u } }$,

$M _ { 1 } ^ { h } = \frac{ M _ { 1 } ^ { y } - η_ { 1 } M _ { 0 } ^ { y } } { M _ { 0 } ^ { u } }$,

$M _ { 2 } ^ { h } = \frac{ M _ { 2 } ^ { y } -2 η_ { 1 } M _ { 1 } ^ { y } + \left \{ 2(η _ { 1 } ) ^ { 2 } -η _ { 2 } \right \} M _ { 0 } ^ { y } } { M _ { 0 } ^ { u } }$,

$M _ { 3 } ^ { h } = \frac{ M _ { 3 } ^ { y } -3 η_ { 1 } M _ { 2 } ^ { y } -3(η_ { 2 } -2 η _ { 1 } ^ { 2 } )M _ { 1 } ^ { y } - \left \{ η_ { 3 } +6(η _ { 1 } ^ { 3 } - η_ { 1 } \Phi _ { 2 } ) \right \} M _ { 0 } ^ { y } } { M _ { 0 } ^ { u } } ,$

$M _ { 4 } ^ { h } = \frac{ M _ { 4 } ^ { y } -4 \left \{ M _ { 1 } ^ { u } M _ { 3 } ^ { h } + η_ { 3 } (M _ { 1 } ^ { y } -η _ { 1 } M _ { 0 } ^ { y } ) \right \} -6M _ { 2 } ^ { u } M _ { 2 } ^ { h } -η_ { 4 } M _ { 0 } ^ { y } } { M _ { 0 } ^ { u } } ,$

여기서, $η _ { k } = \frac{ M _ { k } ^ { u } } { M _ { 0 } ^ { u } } ,~~ k=1,2,3...$.

입력과 출력의 샘플데이터를 취득하면 입출력 모멘트는 (2)에 의해 구할 수 있으므로 (15)로부터 우리는 미지 프로세스의 k차 모멘트를 계산할 수 있다. 따라서 파라미터 식별 공식 유도 시 프로세스의 모멘트 $M _ { k } ^ { h }$를 k차 까지 알고 있다고 간주한다.

3.1 일반적인 이산시간저차모델의 식별

저차 근사화 모델 $ \widehat {H} _L (z ^ { -1 } )$의 모멘트를 $ M _ {k} ^ { \widehat {h} _{ L }} $이라 하면 모멘트 정합 방법은 다음을 만족하도록 $ \widehat { H } _L (z ^ { -1 } )$을 구하는 것이다.

(16)

$M _ {k} ^ { \widehat {h} _{ L }} = M _ { k } ^ { h } , k=0, 1, 2,\cdots, N$

${ \widehat H } _L (z ^ { -1 } )$이 $H(z ^ { -1 } )$의 근사화 모델이면 다음 식으로 나타낼 수 있다.

(17)

$H(z ^ { -1 } ) \cong { \widehat { H } } _ { L } (z ^ { -1 } )= \frac{ N _ { L } (z ^ { -1 } ) } { D _ { L } (z ^ { -1 } ) } z ^ { - { \widehat { d } } }$

$⇔ H(z ^ { -1 } )D _ { L } (z ^ { -1 } ) \cong N _ { L } (z ^ { -1 } ) z ^ { - { \widehat { d } } }$

식 (6),(7)로부터 모멘트 정합조건 (16)은 두 모델의 k차까지 도함수(derivative)를 정합시키는 것과 등가관계이다. 이제 (17)의 두 번째 식을 $z=1$에서 0차부터 4차까지 도함수 식을 전개하기로 한다.

(18)

$H(z ^ { -1 } )D _ { L } (z ^ { -1 } )| _ { z=1 } =N _ { L } (z ^ { -1 } ) z ^ { - { \widehat { d } } } | _ { z=1 } ~$

$⇒ ~H(1)D _ { L } (1)=N _ { L } (1)$

식 (6)을 (18)에 대입하면

(19)

$M_0^h (1+a _ { 1 } +a _ { 2 } + \cdots +a _ { l } ) =(b _ { 1 } +b _ { 2 } +b _ { 3 } + \cdots b _ { q } )$

$z=1$에서 (17)의 1차도함수를 취하고 (6)을 대입하면

(20)

$H ' (1)D _ { L } (1) + H(1)D_ { L } ' (1)=N _ { L } ' (1) + { \widehat d } N _ { L } (1)$

$\Rightarrow ~M _ { 1 } ^ { h } (1+a _ { 1 } +a _ { 2 } + \cdots +a _ { l } )+M _ { 0 } ^ { h } (a _ { 1 } +2a _ { 2 } + \cdots +l a _ { l } )$ $=(b _ { 1 } +2b _ { 2 } + \cdots +q b _ { q } )+ { \widehat { d } } (b _ { 1 } +b _ { 2 } + \cdots +b _ { q } )$

이와 유사하게 $z=1$에서 (17)의 2~4차 도함수를 취하여 전개하면

2차도함수:

(22)

$H '' (1)D _ { L } (1)+ 2H '(1)D_ { L } ' (1) + H(1)D_ { L } '' (1)$

$=N _ { L } '' (1) +2 N _ { L } ' (1) { \widehat d } + N _ { L } (1) { \widehat d } ( { \widehat d } -1)$

$\Rightarrow (M_2^h -M_1^h ) (1+a _ { 1 } +a _ { 2 } + \cdots +a _ { l } )$

$+2 M_1^h (a _ { 1 } +2a _ { 2 } + \cdots +l a _ { l } ) +M_0^h (2a_2 +6a_3 +cdots+l (l-1)a_l ) $

$=(2b_2 +6b_3 +cdots+q (q-1)b_q )+2 (b _ { 1 } +2b _ { 2 } + \cdots +q b _ { q } ) { \widehat d }$

$+ (b _ { 1 } +b _ { 2 } + \cdots+ b _ { q } ) { \widehat d } ( { \widehat d } -1)$

3차도함수:

(23)

$H ''\:' (1)D _ { L } (1)+3H ''(1)D_ { L } ' (1) + 3H '(1)D_ { L } '' (1) + H(1)D_ { L } ''\:' (1) $

$=N _ { L } ''\:' (1) +3 N _ { L } '' (1) { \widehat d } +3 N _ { L } ' (1) { \widehat d } ( { \widehat d } -1)+ N _ { L } (1) { \widehat d } ( { \widehat d } -1)( { \widehat d } -2)$

$\Rightarrow ~(M _ { 3 } ^ { h } -3M _ { 2 } ^ { h } +2M _ { 1 } ^ { h } )(1+a _ { 1 } + \cdots +a _ { l } )+3(M _ { 2 } ^ { h } -M _ { 1 } ^ { h } )(a _ { 1 } + \cdots +l a _ { l } )$

$+3M _ { 1 } ^ { h } (2a _ { 2 } + \cdots +l (l-1)a _ { l } )+M _ { 0 } ^ { h } (6a _ { 3 } + \cdots +l(l-1)(l-2)a _ { l } )$

$=(6b _ { 3 } + \cdots +q(q-1)(q-2)b _ { q } )+3(2b _ { 2 } + \cdots +q (q-1)b _ { q } ) { \widehat { d } }$

$+3(b _ { 1 } + \cdots +q b _ { q } ) { \widehat { d } } ( { \widehat { d } } -1)+ (b _ { 1 } + \cdots +b _ { q } ) { \widehat { d } } ( { \widehat { d } } -1)( { \widehat { d } } -2)$

4차도함수:

(23)

$H ^ { (4) } D _ { L } (1)+4H ''\:' (1)D _ { L } ' (1)+6H '' (1)D _ { L } '' (1)+4H ' (1)D _ { L } ''\:' (1)$

$+H(1)D _ { L } ^ { (4) } (1)$

$=N _ { L } ^ { (4) } (1) +4N _ { L } ''\:' (1) { \widehat { d } } +6N _ { L } '' (1) { \widehat { d } } ( { \widehat { d } } -1)+4 N _ { L } ' (1) { \widehat { d } } ( { \widehat { d } } -1)( { \widehat { d } } -2)$

$+N _ { L } (1) { \widehat { d } } ( { \widehat { d } } -1)( { \widehat { d } } -2)( { \widehat { d } } -3)$

$\Rightarrow ~(M _ { 4 } ^ { h } -6M _ { 3 } ^ { h } +11 M _ { 2 } ^ { h } -6M _ { 1 } ^ { h } )(1+a _ { 1 } + \cdots +a _ { l } )+4(M _ { 3 } ^ { h } -3 M_ { 2 } ^ { h } +2M _ { 1 } ^ { h } )$

$(a _ { 1 } + \cdots +l a _ { l } )+6(M _ { 2 } ^ { h } -M _ { 1 } ^ { h } ) (2a _ { 2 } + \cdots +l (l-1)a _ { l } )+4M _ { 1 } ^ { h } $

$(6a _ { 3 } + \cdots +l(l-1)(l-2)a _ { l } ) +M_0^h (24a_4 +cdots +l(l-1)(l-2)(l-3)a_L )$

$=(24b_4 +cdots +q(q-1)(q-2)(q-3)b_q )+4(6b _ { 3 } + \cdots +q(q-1)(q-2)b _ { q } ) { \widehat { d } }$

$+4(2b _ { 2 } + \cdots +q (q-1)b _ { q } ) { \widehat { d } } ( { \widehat { d } } -1)+4(b _ { 1 } + \cdots +q b _ { q } ) { \widehat { d } } ( { \widehat { d } } -1)( { \widehat { d } } -2)$ $+ (b _ { 1 } + \cdots +b _ { q } ) { \widehat { d } } ( { \widehat { d } } -1)( { \widehat { d } } -2)( { \widehat { d } } -3)$

지연시간 $\widehat { d }$ 가 주어졌다고 가정하고 (19)~ (23)을 미지 파라미터에 대해 행렬 형태로 정리하면

(24)

$\Phi_L \theta=v_L$,

여기서

(25)

$\theta :=[ { \widehat a_1 } { \widehat a_2 } \cdots { \widehat a_l } { \widehat b_1 } { \widehat b_2 } \cdots { \widehat b_q ]^T }$

(26)

$\phi_L = \begin{bmatrix} M _ { 0 } ^ { h } M _ { 0 } ^ { h } \cdots M _ { 0 } ^ { h } -1 -1 \cdots -1 \\ M _ { 1 } ^ { h } +M _ { 0 } ^ { h } M _ { 1 } ^ { h } +2M _ { 0 } ^ { h } \cdots M _ { 1 } ^ { h } +l M _ { 0 } ^ { h } - (\widehat { d } +1) -( \widehat { d } +2) \cdots -( \widehat { d } +q) \\ \phi_ { 31 } \phi_ { 32 } \cdots \phi_ { 3l } - \widehat { d } (\widehat { d } +1) -( \widehat { d } +1) (\widehat { d } +2) \cdots -( \widehat { d } +q-1)( \widehat { d } +q) \\ \phi_ { 41 } \phi_ { 42 } \cdots \phi_ { 4l } \phi_ { 4(l+1) } \phi_ { 4(l+2) } \cdots \phi_ { 4(l+q) } \\ \phi_ { 51 } \phi_ { 52 } \cdots \phi_ { 5l } \phi_ { 5(l+1) } \phi_ { 5(l+2) } \cdots \phi_ { 5(l+q) } \\ \cdots \cdots \end{bmatrix}$

(27)

$v_L= \begin{bmatrix} -M_0^h \\ -M_1^h \\ -(M_2^h -M_1^h ) \\ -(M_3^h -3M_2^h +2M_1^h ) \\ -( M_4^h -6M_3^h +11M_2^h -6M_1^h ) \\ \vdots \end{bmatrix}$

식 (26)의 행렬 요소는 다음과 같다.

(28)

$\phi_ { 31 } =M_2^h +M_1^h$, $\phi_ { 32 } =M_2^h +3M_1^h +2M_0^h ,$

$\phi_ { 3l } =M_2^h +(2l-1)M_1^h +l(l-1)M_0^h ,$

$\phi_ { 41 } =M_3^h -M_1^h ,$ $\phi_ { 42 } =M_3^h +3M_2^h +2M_1^h ,$

$\phi_ { 4l } =M_3^h +3(l-1)M_2^h +[3l(l-2)+2]M_1^h +l(l-1)(l-2)M_0^h ,$

$\phi_ { 4(l+1) } =-( \widehat { d } -1) { \widehat d } (\widehat { d } +1),$

$\phi_ { 4(l+2) } =- \widehat { d } ( { \widehat d } +1) (\widehat { d } +2),$ $\phi_ { 4(l+q) } =-( \widehat { d } +q -2) ( \widehat { d } +q -1) (\widehat { d } +q),$

$\phi_ { 51 } =M_4^h -2M_3^h -M_2^h +2M_1^h ,$ $\phi_ { 52 } =M_4^h +2M_3^h -M_2^h -2M_1^h ,$

$\phi_ { 5l } =M_4^h +(4l-6)M_3^h +[6l(l-3)+11]M_2^h +[4l(l-1)(l-2)$

$-2l(3l-7)-6]M_1^h +l(l-1)(l-2)(l-3)M_0^h ,$

$\phi_ { 5(l+1) } =-( \widehat { d } -2) ( \widehat { d } -1) { \widehat d } (\widehat { d } +1),$ $\phi_ { 5(l+2) } =- ( \widehat { d } -1) \widehat { d } ( { \widehat d } +1) (\widehat { d } +2),$

$\phi_ { 5(l+q) } =-( \widehat { d } +q -3) ( \widehat { d } +q -2) ( \widehat { d } +q -1) (\widehat { d } +q).$

식 (24)의 확장행렬 $[\Phi_L \vdots v_L ]$에 행 연산(row operation)을 수행하면 (26), (27)보다 규칙성을 갖는 다음 선형방정식을 유도할 수 있다.

(29)

${ \bar \Phi } _L \theta= { \bar v } _L$,

여기서,

(30)

${ \bar { \Phi } } _ { L } = \begin{bmatrix} M _ { 0 } ^ { h } M _ { 0 } ^ { h } \cdots M _ { 0 } ^ { h } -1 -1 \cdots -1 \\ M _ { 1 } ^ { h } M _ { 1 } ^ { h } +M _ { 0 } ^ { h } \cdots M _ { 1 } ^ { h } +(l-1) M _ { 0 } ^ { h } - { \widehat { d } } -( { \widehat { d } } +1) \cdots -( { \widehat { d } } +q-1) \\ M _ { 2 } ^ { h } M _ { 2 } ^ { h } +2M _ { 1 } ^ { h } +M _ { 0 } ^ { h } \cdots M _ { 2 } ^ { h } +2(l-1)M _ { 1 } ^ { h } +(l-1) ^ { 2 } M _ { 0 } ^ { h } - { \widehat { d } } ^ { 2 } -( { \widehat { d } } +1) ^ { 2 } \cdots -( { \widehat { d } } +q-1) ^ { 2 } \\ \vdots \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \vdots \: \: \: \: \: \: \: \: \: \: \cdots \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \vdots \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \vdots \: \: \: \: \: \: \: \: \: \: \: \: \: \: \: \vdots \: \: \: \: \: \: \: \: \: \: \cdots \: \: \: \: \: \: \: \: \: \: \vdots \\M _ { N } ^ { h } \: \: \: \: \: \: \: \: \: \: { \bar { \phi } } _ { (N+1)2 } \: \: \: \: \: \: \: \: \: \: \cdots \: \: \: \: \: \: \: \: \: \: { \bar { \phi } } _ { (N+1)l } \: \: \: \: \: \: \: \: \: \: \cdots \: \: \: \: \: \: \: \: \: \: \cdots \: \: \: \: \: \: \: \: \: \: \cdots \: \: \: \: \: \: \: \: \: \: { \bar { \phi } } _ { (N+1)(l+q) } \end{bmatrix}$

식 (30)의 마지막 행의 요소는

(31)

${ \bar { \Phi } } _ { (N+1)2 } = \displaystyle \sum _ { k=0 } ^ { N } \begin{pmatrix}N\\i\end{pmatrix} M_k^h ,$

${ \bar { \Phi } } _ { (N+1)l } = \displaystyle \sum _ { k=0 } ^ { N } \begin{pmatrix}N\\i\end{pmatrix} (l-1)^ { N-k } M_k^h , $

${ \bar { \Phi } } _ { (N+1)(l+q) } = -( { \widehat d } +q-1)^N .$

(32)

${ \bar v } _L= \begin{bmatrix} -M_0^h \\ -(M_1^h -M_0^h ) \\ -(M_2^h -2M_1^h +M_0^h )) \\ -(M_3^h -3M_2^h +3M_1^h -M_0^h ) \\ \vdots \\ - \sum _ { N } ^ { k=0 } (-1)^ { k+1 } \begin{pmatrix}N\\i\end{pmatrix} M_k^h \end{bmatrix}$

미지 파라미터 수를 $p=l+q$ 라 하면 (23)처럼 $N=p-1$차까지 모멘트 정합이 이루어지도록 수식을 전개하면 (30)의 행렬은 정방행렬이 되고, 각 열벡터가 1차 독립이면 정칙행렬이다. (30)~ (32)로부터 행렬 ${ \bar Phi } _L$, ${ \bar v } _L$는 플랜트의 모멘트와 지연시간으로 구성되기 때문에 저차모델 ${ \widehat H } _L (z ^ { -1 } )$의 파라미터는 (29)의 유일해로 구해진다.

부언 1: 식 (30)의 ${ \bar { \Phi } }_L$이 정칙이기만 하면 (29)는 해가 존재한다. 이 조건은 $M_0^h$, $M_1^h, M_2^h, \cdots ,M_N^h$가 서로 다르기만 하면 $\widehat d$값에 관계없이 성립함을 알 수 있다. 식 (15)에 의해 입출력 모멘트로부터 미지 프로세스의 모멘트가 정확하게 구해진다면 시험입력에 지속여기 조건은 필요하지 않다.

3.2 FOTD, SOTD, SODF 모델의 식별 공식

이제 (29)~(32)의 결과로부터 세 가지 저차 모델, FOTD, SOTD, SODF의 파라미터를 식별하는 공식을 나타내기로 한다. 식 (12)로부터 FOTD 모델:

$FOTD : { \widehat H } _ { FD } (z ^ { -1 } )= \frac{ \widehat { b _ { 1 } } z ^ { -1 } + \widehat { b _ { 2 } } z ^ { -2 } } { 1+ \widehat { a _ { 1 } } z ^ { -1 } } z ^ { - \widehat { d } }$

에 대해 파라미터 벡터를 $ \theta _ { FD } := [ \widehat a _1 \widehat b _1 \widehat b _2 ]^T $라 정의하면, 미지 이산시간 프로세스에 대해 모멘트정합 조건을 만족하는 지연시간 $\widehat { d }$인 FOTD 모델의 파라미터는 다음 식으로 결정된다.

(33)

$ \theta_ { FD } = { \bar \Phi } _FD^ { -1 } { \bar v } _ { FD } ,$

여기서

(34)

${ \bar \Phi } _ { FD } = \begin{bmatrix} M _ { 0 } ^ { h } \: \: \: \: \: -1 \: \: \: \: \: -1 \\ M _ { 1 } ^ { h } -\widehat { d } -( \widehat { d } +1) \\ M _ { 2 } ^ { h } - { \widehat d } ^2 -( \widehat { d } +1)^2 \end{bmatrix} $,

(35)

${ \bar v } _ { FD } = \begin{bmatrix} -M_0^h \\ -(M_1^h -M_0^h ) \\ -(M_2^h -2M_1^h +M_0^h ) \end{bmatrix} .$

FOTD 모델은 3개의 파라미터를 식별하므로 프로세스의 모멘트는 $N=p-1=(l+q)-1=2$차 까지 구하면 된다.

마찬가지로 SOTD 모델:

$SOTD : \widehat H _ { SD } ( z ^ { -1 } )= \frac { \widehat { b } _ { 1 } z ^ { -1 } + \widehat { b } _ { 2 } z ^ { -2 } + \widehat { b } _ { 3 } z ^ { -3 } } { 1+ \widehat { a } _ { 1 } z ^ { -1 } + \widehat { a } _ { 2 } z ^ { -2 } } z ^ { - \widehat { d } }$

에 대해 파라미터 벡터를 $\theta_ { SD } :=[ { \widehat a_1 } { \widehat a_2 } { \widehat b_1 } { \widehat b_2 } { \widehat b_3 } ]^T }$라 정의하면, 미지 이산시간 프로세스에 대해 모멘트정합 조건을 만족하는 지연시간 $\widehat { d }$인 SOTD 모델의 파라미터는 다음 식으로 결정된다.

(36)

$theta_ { SD } = { \bar Phi } _ { SD } ^ { -1 } { \bar v } _ { SD } ,$

여기서

(37)

${ \bar \Phi } _ { SD } = \begin{bmatrix} M _ { 0 } ^ { h } M _ { 0 } ^ { h } -1 ~-1 ~-1 \\ M _ { 1 } ^ { h } M _ { 1 } ^ { h } +M _ { 0 } ^ { h } - { \widehat d } ~-( \widehat { d } +1) ~-( \widehat { d } +2) \\ M _ { 2 } ^ { h } M _ { 2 } ^ { h } +2M _ { 1 } ^ { h } +M _ { 0 } ^ { h } - { \widehat d } ^2 ~-( \widehat { d } +1)^2 ~-( \widehat { d } +2)^2 \\ M _ { 3 } ^ { h } M _ { 3 } ^ { h } +3M _ { 2 } ^ { h } +3M _ { 1 } ^ { h } +M _ { 0 } ^ { h } - { \widehat d } ^3 ~-( \widehat { d } +1)^3 ~-( \widehat { d } +2)^3 \\ M _ { 4 } ^ { h } M _ { 4 } ^ { h } +4M _ { 3 } ^ { h } +6M _ { 2 } ^ { h } +4M _ { 1 } ^ { h } +M _ { 0 } ^ { h } - { \widehat d } ^4 ~-( \widehat { d } +1)^4 ~-( \widehat { d } +2)^4 \end{bmatrix} ,$

(38)

${ \bar v } _ { SD } = \begin{bmatrix} -M_0^h \\ -(M_1^h -M_0^h ) \\ -(M_2^h -2M_1^h +M_0^h ) \\ -(M_3^h -3M_2^h +3M_1^h -M_0^h ) \\ -(M_4^h -4M_3^h +6M_2^h -4M_1^h +M_0^h ) \end{bmatrix} .$

SOTD 모델은 5개의 파라미터를 식별하므로 프로세스의 모멘트는 $N=(l+q)-1=4$차까지 구하면 된다.

SODF모델은 (37)에서 $\widehat { d } =0$로 놓으면 (36)으로부터 바로 구해진다.

식 (14) 모델에서 분모와 분자의 차수를 임의로 선택해도 (30)과 (32)에서 각 파라미터에 대응하는 행과 열을 남겨두고 나머지는 제거하는 방식으로 (29)를 수정하면 다른 형식의 저차모델에도 동일하게 적용할 수 있다.

부언 2: 연속시간 시스템의 모멘트 정합에 의한 연속시간 저차모델링 방법[7]으로 연속시간 FOTD, SOTD 모델을 구하여 이산화하면 (33), (36)으로 직접 구한 모델과 거의 유사한 결과를 낳는다. 그런데 연속시간 시스템의 모멘트를 정확하게 구하려면 입・출력 데이터를 가능한 많이 필요로 하는 반면 이산시간 시스템의 모멘트는 샘플 시간에서의 데이터만 이용하기 때문에 계산이 간단하고 정확도도 높다.

3.3 지연시간의 추정

이산시간 전달함수를 모델링할 때 지연시간 $\widehat { d }$는 정수이어야 한다. 여기서 제시하는 지연시간 추정방법은 먼저 지연시간 후보 집합을 선택하고 이 집합의 각 지연시간에 대응하는 저차모델 파라미터를 3.1, 3.2절에서 제시한 공식을 이용하여 구한다. 그리고 이 모델 중에서 실제 펄스응답과의 IAE를 비교하여 최소 IAE를 주는 지연시간과 파라미터를 최종 식별모델로 결정하게 된다.

펄스응답의 초기 응답 지연시간을 $L_o$[sec] $L_o =d_j T_s + \lambda T_s , ~(0 \leq \lambda \leq 1)$에 의해 $d_j$를 구한다. 이 $d_j$를 중심 값으로 $±T_s$, $±2T_s$인 정수 값을 다음과 같이 지연시간 후보 집합으로 정의한다.

(40)

$D:= \left \{ d | \left [ d _ { j-2 } , d _ { j-1 } , d _ { j } , d _ { j+1 } , d _ { j+2 } \right ] \right \}$

여기서 $d_ { j+-k } =d_j +-k T_s ,~k=0, 1, 2.$ 미지 플랜트 $H(z ^ { -1 } )$의 펄스응답, 스텝응답, 임펄스 응답을 각각 ${ y _ { p } } (kT)$, ${ y } _ { i } (kT)$, ${ y } _ { s } (kT)$라 하고, 식별한 모델 ${ \widehat H } _L (z ^ { -1 } )$의 펄스응답, 스텝응답, 임펄스 응답을 각각 $\widehat { y _ { p } } (kT)$, $\widehat { y } _ { i } (kT)$, $\widehat { y } _ { s } (kT)$라 표기하기로 한다. 각 경우의 IAE는 다음과 같이 정의한다.

(40)

$J _ { p } := \sum _ { k=0 } ^ { \infty } |y _ { p } (kT _ { s } )- \widehat { y _ { p } } (kT _ { s } )|,$

이제 펄스응답의 IAE를 최소화시키는 지연시간 추정은 다음 식으로 나타낼 수 있다.

(41)

${ \widehat { d } } =arg \: min(J _ { p } (d)) , { \forall } d \in D$

위 주요결과와 2장의 예비결과를 요약하면 다음 알고리즘으로 정리된다.

알고리즘 1:

[스텝로부터 를 구하고 (39)의 지연시간 후보 집합을 정한다. 식 (14)에서 식별할 저차 모델의 분모 차수($l$), 분자 차수($q$)를 선택한다. $N=l+q-1$이다. [스텝$N$차까지 프로세스의 모멘트$M_0^h$, $M_1^h$, $M_2^h \cdots, M_N^h$를 계산한다. [스텝에 대해 (29)(FOTD, SOTD이면 (33), (36))를 이용하여 각 지연시간에 대응하는 저차모델을 구한다. [스텝 [스텝

[스텝1] 프로세스에 단일 펄스(구형파, 반주기 정현파 등) 입력을 인가하여 입력과 출력의 샘플데이터를 취득한다. [샘플링 시간 $T_s$ ]

[스텝2] 펄스응답의 초기 응답 지연시간 $L_o$로부터 $d_j$를 구하고 (39)의 지연시간 후보 집합을 정한다. 식 (14)에서 식별할 저차 모델의 분모 차수($l$), 분자 차수($q$)를 선택한다. $N=l+q-1$이다.

[스텝3] 식 (2)를 이용 입출력의 모멘트를 계산하고, (15)를 이용하여 $N$차까지 프로세스의 모멘트 $M_0^h$, $M_1^h$, $M_2^h \cdots, M_N^h$ 를 계산한다.

[스텝4] ${ \widehat d } ~in D$에 대해 (29) (FOTD, SOTD이면 (33), (36))를 이용하여 각 지연시간에 대응하는 저차모델을 구한다.

[스텝5] [스텝 4]에서 구한 각 모델에 대해 시뮬레이션으로 시험입력과 동일한 펄스입력을 인가하여 응답을 구하고 (40)의 IAE를 계산한다. (41)로부터 최소 IAE를 낳는 지연시간과 대응하는 모델을 최적 저차모델로 결정한다.

[스텝6] 식별할 모델의 차수를 너무 낮게 선택하면 응답일치성이 떨어질 수 있음. [스텝 2]로 돌아가 저차 모델의 분모 차수($l$), 분자 차수($q$)를 다시 선택하고 위 과정을 반복한다.

4. 시뮬레이션

이 절에서는 한 예제를 통해 어떻게 알고리즘 1을 적용하여 저차 모델을 식별하는가를 보여주고자 한다.

예제 1: 다음 전달함수를 미지 프로세스로 고려한다.

$G(s)= \frac{ 1 } { (2s+1)(s+1)(0.5s+1)(0.1s+1) }$

샘플링 시간 $T_s =0.2 sec$일 때 ZOH를 포함한 이산시간 플랜트 모델은

(42)

$H(z ^ { -1 } )= \frac{ 0.405 { \rm x } 10^-3 z ^ { -1 } +2.786 { \rm x } 10^-3 z ^ { -2 } +1.646 { \rm x } 10^-3 z ^ { -3 } +8.087 { \rm x } 10^-5 z ^ { -4 } } { 1-2.5292z ^ { -1 } +2.2201z ^ { -2 } -0.7532z ^ { -3 } +0.0672z ^ { -4 } }$

식 (42)의 4차 이산시간 프로세스를 펄스응답으로부터 알고리즘 1을 적용하여 SOTD 모델을 구해보기로 한다. < 그림 2>는 프로세스 (42)에 인가한 구형파 펄스 입력과 출력이다. < 그림 2 (b)>로부터 초기 응답지연 시간이 약 $L_o =0.4$초로 나타났다. $d _ { j } =0.4/T_S =2$가 되어 가능한 지연시간 집합을 $D = \left \{ { \widehat d } | 0,1,2,3,4 \right \}$로 정한다. 식 (13)의 SOTD모델은 분모 분자 차수가 각각 $l=2$, $q=3$이므로 $N=(l+q)-1=4$ 차까지 프로세스의 모멘트를 구하면 된다.

그림2 구형파 펄스입력과 응답: (a) 구형파 펄스 입력 $u(k)$, (b) 펄스 응답 $y(k)$ .

Fig.2 Rectangular pulse input $u(k)$ and response $y(k)$.

../../Resources/kiee/KIEE.2018.67.8.1062/fig2.png

< 표 1>은 입출력 모멘트와 (15)를 이용하여 구한 플랜트의 모멘트 값이다.

지연시간을 추정하기 위해 ${ \widehat d } ~ \in D$ 에 대해 (36)~ (38)을 이용하여 5개의 모델을 구하고 IAE를 계산하면 < 표 2>와 같다. $\widehat { d } =2$ 일 때 IAE 값이 최소가 되고 (42)에 대해 다음 최적 SOTD모델을 얻었다.

(43)

${ \widehat H } _ { SD } (z ^ { -1 } )= \frac{ 0.03091z ^ { -1 } -0.04829z ^ { -2 } +0.03369z ^ { -3 } } { 1-1.7338z ^ { -1 } +0.7501z ^ { -2 } } z ^ { -2 }$

< 그림 3>은 (42)의 플랜트$H(z ^ { -1 } )$와 SOTD 식별 모델 ${ \widehat H } _ { SD } (z ^ { -1 } )$의 시간응답을 비교한 것이다. 4차를 2차로 모델링하였음에도 두 모델의 임펄스 응답과 스텝응답이 매우 유사한 응답을 보여주고 있다. < 그림 3(a)>의 임펄스 응답초기에 두 응답 간 차이가 나타나는 것은 4차인 원래 시스템을 2차로 모델링한데서 기인한다. 이는 (42)를 직접 가장 정밀도 높은 차수 줄임 기법을 적용하여 2차로 줄인 모델에서도 유사한 응답 차이를 쉽게 확인할 수 있다. 이 예는 제안한 방법에 따라 (42)의 4차 플랜트에 매우 간단한 구형파 펄스 입력을 시험입력으로 인가하고 그 응답으로부터 (43) 모델의 5개 파라미터를 식별하였음에도 잘 식별할 수 있음을 보여주고 있다.

그림3 예제 1의 4차 플랜트 $H( z^-1)$와 SOTD 식별 모델 $ \widehat H _{ SD }( z^-1)$의 시간응답 비교

Fig. 3 Impulse and step responses of plant $H( z^-1)$ and the identified SOTD model $ \widehat H _{ SD }( z^-1)$ in Example 1.

../../Resources/kiee/KIEE.2018.67.8.1062/fig3.png

< 표 1>에서 구한 플랜트 모멘트 값을 (33)과 알고리즘에 적용하면 바로 FOTD 모델이 얻어진다. 여기서는 생략한다.

표 1 입・출력과 시스템의 모멘트 값

Table 1 Values of input, output, and system moments

모멘트 차수

입력모멘트 ( $ M_k^u $ )

출력모멘트

( $ M_k^y $ )

플랜트모멘트 ( $ M_k^h $ )

0

5.0000

5.0000

1.0000

1

10.0000

102.5000

18.5000

2

30.0000

2769.1667

473.8333

3

100.0000

96000.000

15916.0000

4

354.0000

4090182.4864

672099.6973

표 2 ${\widehat{d}} \in D$ 에 대한 SOTD모델의 IAE 값

Table 2 IAE values of the SOTD models for ${\widehat{d}} \in D$

지연시간

( $ \widehat d $ )

0

1

2

3

4

IAE

0.0414

0.0273

0.0186

0.0235

0.0561

5. 결 론

본 논문에서는 간단한 펄스응답으로부터 직접 이산시간 저차모델을 식별하는 알고리즘을 제시하였다. 기본 개념은 플랜트의 모멘트와 식별한 저차 모델의 모멘트가 0차부터 $N$차까지 서로 같아지도록 모델의 파라미터를 결정하는 것이다. 제시한 방법은 모멘트 정합에 의한 차수 줄임과 모델 파라미터 식별을 조합한 것이라 할 수 있다. 이를 위해 유한 크기의 입력과 그 응답으로부터 미지 선형 시스템의 모멘트를 구할 수 있음을 해석적으로 보였다. 지연시간 추정과 모델 파라미터를 분리하여 식별하였다. 먼저 지연시간이 주어진 조건에서 모멘트 정합에 근거한 일반적인 저차모델의 식별 공식을 유도하였으며 이 결과로부터 실제 모델링 시 가장 많이 나타나는 FOTD, SOTD 모델을 식별하는 공식을 제시하였다. 저차모델의 지연시간은 가능한 지연시간 집합을 설정하고 각 지연시간에 대응하는 모델을 구하고 이 모델 집합으로부터 펄스응답을 구하여 IAE를 최소화 시키는 지연시간을 결정하도록 하였다.

식별공식의 해의 존재 조건으로부터 입력에 지속여기 조건을 요구하지 않음을 보였다. 시험입력은 구형파 펄스, 반주기 정현파와 같은 간단한 형태라도 충분하기 때문에 실제 문제에 적용하기가 매우 용이하다.

예에서 보여주듯이 모멘트는 고차로 갈수록 그 값이 커지기 때문에 식별할 파라미터 수가 많아지면 행렬 ${ \bar \Phi } _L$의 역행렬을 구할 때 수치적 문제가 생길 수 있다. 따름정리 1 (또는 식 (15))에 따라 계산하는 시스템의 모멘트의 정밀도가 시험입력의 파형, 출력에 포함되는 잡음 등에 의해 어떤 영향을 받는지는 추후 과제로 남겨둔다.

감사의 글

이 논문은 2015년도 정부(교육부)의 재원으로 한국연구재단의 지원을 받아 수행된 기초연구사업임(No. NRF- 2015R1D1A01060997)

References

1 
K. J. Åström, T. Hägglund, , Advanced PID Control, ISA-Instrumentation, Systems, and Automation Society, 2006
2 
L. Ljung, , System Identification: Theory for the User, Prentice-Hall, Inc., Upper Saddle River, New Jersey, 1987
3 
T. Söderström, P. Stoica, , System Identification ,Prentice-Hall International (UK) Ltd, 1989
4 
L. Pernebo, L. M. Silverman, "Model Reduction via Balanced State Space Representations", IEEE Transactions on Automatic Control, Vol. 27, No. 2, pp. 382-387, 1982DOI
5 
A. C. Antoulas, , Approximation of Large-scale Dynamical Systems, SIAM, Philadelphia, 2005
6 
A. Astolfi, "Model reduction by moment matching for linear and nonlinear systems", IEEE Transactions on Automatic Control, Vol. 55, No. 10, pp. 2321-2336, 2010DOI
7 
Y. C. Kim, L. Jin, "Robust identification of continuous-time low-order models using moments of a single rectangular pulse response", Journal of Process Control, Vol. 23, No. 5, pp. 682-695, 2013DOI
8 
J. E. Lee, Y. C. Kim, "Experimental modeling of acceleration and brake systems for autonomous vehicle", Transactions of the Korean Institute of Electrical Engineers, Vol. 65, No. 4, pp. 642-651, 2016DOI

저자소개

황 지 호 (Jiho Hwang)
../../Resources/kiee/KIEE.2018.67.8.1062/au1.png

2016. 2 충북대학교 전자공학부 공학사

2018.2 충북대학교 대학원 전자공학전공 석사

관심분야: 시스템 모델링 및 디지털 제어 시스템 설계

차 승 표(Seungpyo Cha)
../../Resources/kiee/KIEE.2018.67.8.1062/au2.png

2017. 2 충북대학교 전기공학부 공학사

2018. 현재 충북대학교 대학원 전자공학전공 석사과정 재학 중

관심분야: 시스템 모델링 및 디지털 제어 시스템 설계

김 영 철(Young Chol Kim)
../../Resources/kiee/KIEE.2018.67.8.1062/au3.png

1981.2 고려대학교 전기공학과(공학사)

1983.2 서울대학교 전기공학과(공학석사)

1987.8 서울대학교 전기공학과(공학박사)

1988년~현재 충북대학교 전자공학부 교수 겸 컴퓨터정보통신연구소 연구원

1992.1 ~1993.1 미국 Texas A & M Univ. Post-Doctoral Fellow

2001.1~2002.1 미국 Vanderbilt Univ./ Tennessee State Univ. 방문교수

2004.1~2008.12 KIEE 제어계측연구회장

2009.1~ 2010.12 KIEE 정보 및 제어 부문회 회장, KIEE 부회장

1988-현재 충북대

관심분야: 저차제어기 설계, 자율 주행차 제어시스템, 시스템 모델링 등