2.1 PMSM 기본식 및 비선형시스템 해석

PMSM의 동작에 관련된 기본 식과 등가회로는 다음과 같이 나타낼 수 있다.

그림. 2. PMSM의 등가회로. (a) d축 등가회로, (b) q축 등가회로.

Fig. 2. The equivalent circuit of PMSM, (a) d-axis equivalent circuit, (b) q-axis equivalent circuit.

여기서, 1) 상태 변수는 $i_{d}$, $i_{q}$ 및 $\omega$는 각각 d-q 축 고정자 전류 및 전동기 각속도를 나타낸다. 반면, 2) 시스템 입력 변수는 $u_{d}$, $u_{q}$이고, 이는 d-𝑞 축 고정자 전압 성분이다. 그리고 $T_{m}$는 외부 부하 토크이다. 파라미터 3) $L_{d}$, $L_{q}$ 및 $R_{s}$는 d-𝑞 축 고정자 인덕턴스 및 고정자 권선 저항이다. 기계적, 자기적 구성요소 4) $\lambda_{af}$, $D$, $J$ 및 $n_{p}$는 각각 영구자석 자속, 점성 감쇠 계수, 극 관성 모멘트, 극쌍 수를 나타낸다. 그림 2는 각각 PMSM 견인전동기의 d축과 q축 등가회로이며, 기본식 식 (1)과 식 (2)의 식으로 나타낼 수 있다. 식 (3)은 전동기의 기본식인 $J\dot\omega =\dfrac{n_{p}P_{e}}{\omega}-T_{m}-T_{D}$를 전개한 형태로, 전기적 입력 $P_{e}$는 그림 2전동기 회로의 구성에 따라 전개할 수 있으며, 감쇠토크는 $T_{D}= D\omega$와 같이 나타낼 수 있다. 위 식은 시스템을 결정하는 내외부 파라미터가 복잡하게 연결되어 있으므로, 시스템(부품)의 내재한 특성을 알기 어렵다. 그래서 ⁽³⁾에서는 아래 식과 같이 무차원화 (Non-dimensionalization)한 식을 제시하며, 원 식의 모든 성질은 보존된다. 이같이 식을 무차원화 (시간 스케일링 기법 및 기하학적 요소 보존기법 등 활용하여 정리)를 함으로써, 철도차량 PMSM의 내부적인 성질을 수학적으로 확인할 수 있다는 것이다.

여기서 윗첨자 *는 식 (1) ~ (3)의 입력변수 및 상태변수를 무차원화를 통해 변환한 변수를 나타낸다.

2.2 분기 해석 (Bifurcation analysis)

분기 해석이란 큰 복잡성을 가진 비선형시스템을 해석하는 일반적인 동역학적인 해석 방법이다. 간단히 말해 분기란, 다음과 같이 미분방정식 형태로 나타낼 수 있는 시스템이 주어진 경우,

여기서, $x$는 상태변수를, $\lambda$는 시스템 파라미터를 나타낸다.초기치 또는 시스템 내 파라미터가 변화할 때, 해당 시스템이 질적으로 변화 (해의 성질) 하는 경우로 정의되고 있다. 다음은 XPPAUT로 확인한 철도 차량용 전동기의 분기 해석 결과이다. 여기서 해당 문헌에서는 $\gamma$를 파라미터로 설정하고, 이를 변화함에 따라서 분기해석을 수행하는데, 이 파라미터는 $\gamma = n_{p}\lambda_{af}^{2}/R_{s}D$ 로 나타낼 수 있으며, 전동기의 자속의 제곱에 비례하는 값이다. 즉, 자속이 커질수록 선형적으로 전류가 세지는 관계가 아닌, 비선형적인 특성을 갖는다는 것을 나타낸다.

그림. 3. 입력 전압 및 부하가 없을 때, PMSM의 평형점 해석 결과.

Fig. 3. Equilibrium analysis of PMSM for zero voltage, zero external torque.

그림 3는 무차원화한 PMSM으로부터 평형점 (Equilibrium)을 찾았을 때, 파라미터의 변화에 따라서 고유치 (Eigenvalue) 해석을 통해 안정도 해석을 한 그림이다. 이 때, 평형점은 각 상태변수가 시간에 따라 변화하지 않을 때의 해로 정의할 수 있으며, 물리적으로는 $i_{d}^{*}$, $i_{q}^{*}$ 및 $\omega^{*}$가 일정한 값을 가질 때를 말한다. 또한, PMSM의 수학적 특성만 도출하기 위하여 수식적으로 토크제어 또는 속도제어가 제외된 형태로, 입력 전압 및 부하토크 또한 고려하지 않았다. 여기서 붉은 색 선은 안정한 해를 갖는 궤적을 나타내며, 검은 색 선은 불안정한 해를 나타낸다. 단 여기서 더 나아가 6, 9 번 포인트는 supercritical Hopf-bifurcation 현상이라 하여, 고유치가 좌반면에서 우반면으로 바뀌는 지점들이다.

그림. 4. 입력 전압 및 부하가 없을 경우, PMSM의 분기현상과 Hopf 분기 지점.

Fig. 4. Hopf bifurcation of PMSM for zero voltage, zero external torque.

그림 4는 Hopf bifurcation point로부터 다시 분기하여 Floquet multiplier analysis (진동, limit cycle)에 대한 해석을 통해 언제 해의 안정한 동요 특성을 잃는지 해석을 수행한 것이다 (파란색 마커). 본 논문에서는 이러한 현상들을 입력 (인가전압, 외부 토크) 이 없을 때를 해석하여 결과를 나타내었지만, 입력이 있는 상태에 서도 동일하게 해석을 수행할 수 있으며, 결과는 참고논문 ⁽³⁾에 포함되어 있다.

요약하자면, 결과의 양상은 크게 다르지 않으나, 특정 수치적 안정도 관련 현상들이 나타나는 지점들이 달라질 수 있다. 다만, 본 논문에서는 추가로 비선형시스템 해석을 직접적으로 수행하지 않는다.

다음 절에서는 실시간 시뮬레이터 구축 모델 내에서 시계열 데이터 분석과 및 통계적인 기법으로 계산하여 시스템 수치적 안정도를 추정하여 검토하는 기술을 논하고자 한다. 먼저 다음 절에서 사례 연구 대상을 소개하도록 한다.

3.1 철도차량 PMSM 제어시스템 실시간 시뮬레이션 성능 검토

그림 5는 OPAL-RT를 활용하기 위해 MATLAB Simulink에서 구성한 FPGA 기반 시뮬레이션 구성도이다. 본 시뮬레이션에서 고려한 사항은 단위 PMSM 1대에 해당하며, OPAL-RT에서 제공한 PMSM 전동기 중 VDQ 모델을 활용하였다. 또한, 시뮬레이션에 활용한 PMSM 전동기의 상세 파라미터는 표 1과 같다.

실시간 시뮬레이션을 수행하기 위해서는 시뮬레이터가 정확히 내부변수 및 출력을 실제와 같은 시간에 출력할 수 있도록 해야 한다. 이는 실시간 시뮬레이션 조건을 달성하기 위해서 시간격이 오프라인 시뮬레이션에 비해 길게 설정되는 이유이다. 실시간 시뮬레이션 단위 시간 간격 간 구성이 표 2에 제공되어 있으며, 그림 6은 이를 간략하게 나타낸 그림이다.

그림. 5. OPAL-RT 시뮬레이터 기반 PMSM과 인버터의 속도제어 구현.

Fig. 5. Configuration of PMSM with inverter and speed controller for OPAL-RT simulation.

시뮬레이션 시간격 $T_{s}$가 주어진 경우, 이 시 간격이 모두 다 연산에 사용되는 것은 아니다. 다시 말하면, 실질적으로 연산 시간 자체는 이 시뮬레이션 시 간격보다 더 짧다. 왜냐하면, 실질적으로 실시간 디지털 시뮬레이터 내부에는 TCP/IP 프로토콜로 코어간 통신 요소가 있으며, HILS의 경우 실제 하드웨어 제어기의 입출력 시간을 고려했을 때도 시뮬레이션의

실시간성이 확보되어야 하기 때문이다. 기본적으로 연산 시간 (Execution cycle)은 선‧후 연산 수행, 처리 시간과 Major 및 minor 연산 시간으로 구성된다. 또한, 실시간 시뮬레이션에서 Idle time 이란 다음 시 간격 순서에 대한 식을 연산하기 전 사용되는 시간으로 정의되며, TCP/IP 통신, 입출력 시간 및 시뮬레이션 클락 대기 시간으로 구성된다. 그러나, 만약 연산 시간을 고려한 시뮬레이터의 동작 시간이 다음 고정 시간 스텝을 초과해 버린다면, 실시간 시뮬레이션은 잘못된 결과를 가져오고, 이를 Overrun이라 말한다 ⁽⁴⁾. 따라서, 시뮬레이션 시 간격 선정은 실시간 시뮬레이션 구현을 위해 어려운 문제이면서도 매우 중요한 결정 요소 중 하나이다.

3.2 PMSM 고유 특성 기반 전동기 부하에 따른 영향 분석

일반적으로 철도차량 실시간 시뮬레이션을 진행할 때, 부하에 대한 요소는 열차 주행 성능 시뮬레이션 (TPS) 기반으로 상세 동적 부하 모델을 사용한다. 실시간 시뮬레이션에서도 물론 정확한 모델링을 전제로 환경 구성을 하지만, 본 논문에서는 추진시스템을 하나의 안정도 추정 ‧ 평가 대상으로서 방법론을 제시하기 위해 간단한 관성 부하를 가정하였으므로, 실시간 시뮬레이션에 따른 실제와의 정합성 등 결과의 수치적 검증을 논하지는 않는다.

3.3 실시간 시뮬레이션 수행 결과 기반 시스템 제어 안정도 추정 기술

위에서 논한 시스템을 활용하여 실시간 시뮬레이션을 수행한

그림. 6. 실시간 디지털 시뮬레이터의 연산시간의 개념도.

Fig. 6. Real time digital simulator time frame.

뒤, 입력 토크를 1) 기울기 2 ($N\bullet m/s$) 램프 및 2) 2 ($N\bullet m$)씩 증가하는 계단 함수를 적용함에 따라 추정된 안정도를 추정하였다. 이는 앞서 2절에서 다룬 이론적인 내용에 해당하며, 그림 7은 외부 토크가 증가함에 따라서 실시간 시뮬레이션 결과로써 시간에 따른 불안정과 상태변수인 $i_{D}$와의 관계 등에 대한 양상을 각각 나타낸다.

그림. 7. 외부 토크($N\bullet m$)를 증가함에 따른 PMSM 실시간 시뮬레이션 상에서 관찰 결과. (a) 기준값 토크 대비 토크의 증가 양상 (램프함수), (b) 램프함수로 설정했을 때의 $i_{d}$ 변화 양상, (c) 기준값 토크 대비 실제 토크의 증가 양상(계단함수) 및 (d) 계단함수로 설정했을 때의 $i_{d}$ 변화 양상.

Fig. 7. Observation of real time simulation PMSM with increasing external torque, (a) Ramp function increment comparing with reference torque, (b) Ramp function increment in $i_{d}$, (c) Stair function increment comparing with reference torque, and (d) Stair function increment in $i_{d}$.

본 논문에서 다루는 사례는 자속기준 제어(FOC ; Field oriented controller) 이므로, $i_{D}$를 0으로 제어하여 최대의 토크를 내는 것이 목적이다. 다만, 그림 7(b)와 그림 7(d)를 살펴보면, 앞서 다루었던 분기 현상을 관찰할 수 있다. 그림 7(b)는 초기에 0으로 제어하는 $i_{D}$가 어느 순간 진동하며 커지며, 최종적으로는 제어계통이 붕괴하는 양상이다. 그림 7(d)의 경우에는 1초마다 값을 변화시키므로, 평형점이 다수가 있다는 것을 알 수 있다.

이렇게 실시간 시뮬레이터에서 취득된 결과는 큰 규모의 시계열 데이터이기 때문에 시계열 데이터만을 기반으로 전반적인 동적시스템의 안정도를 추정할 수 있는 방법의 하나인 시계열 모형 기반 시스템 추정 기법을 도입하여 분석하였다.

3.3.1 시계열 자기회귀 모형 기반 설계모델 분석

자기회귀 (Auto Regressive, AR) 시계열 모형은 전달함수 형태로 나타낼 수 있으며, 즉 시스템 응답을 통해 수치적 안정도를 추정할 수 있다는 것이다. 또한, 관련된 기초 수학 연구들이 이전에 발표된 바 있다 ⁽⁸⁾. 한편, 참고문헌 ⁽⁹⁾에서 주로 시스템의 수치적인 안정성을 검토하기 위하여 AR 모형을 기반으로 한 시스템 고유치 추정을 수행하였다.

3.3.1.1 시계열데이터 정상화 (Stationarize)

PMSM에서 출력되는 상태변수 $i_{D}$, $i_{Q}$, $\omega$의 시계열 분포를 관찰하면, 외부 토크를 지속적으로 증가시키기 때문에 평균이 증가하며, 분산이 증가하는 형상이므로, 데이터의 정상성이 확보되지 않아 정상 시계열 데이터를 전제로 한 자기회귀(Auto regression, AR), 이동평균(moving average, MA) 및 이들의 변형 모형인 ARMA, ARIMA 등을 적용할 수 없다. 따라서 비정상성을 보이는 시계열 데이터를 분석하기 위하여 그림 8과 같이 정상화 알고리즘을 구성하였다. 단, 일반적으로 분산 정상화를 위해 로그만을 취할 때 정상화가 되지만, 짧은 시 간격, 빠른 진동 현상 등으로 인해 정상화가 제대로 이루어지지 않아, 시계열 데이터의 포락선 (envelope)를 취한 뒤, 이에 대한 로그를 취하였음을 밝힌다.

그림. 8. 시계열 정상화 블록 구현.

Fig. 8. Time-series stationarizing block.

그림 9(a)를 보면, $i_{D}$만 분석을 수행하였는데, 이는 FOC 제어의 경우 $i_{D}$를 0으로 하는 제어이기 때문에, 시스템의 안정성을 살펴보기에 충분하다고 판단하였기 때문이다. 따라서 먼저 정상화를 시키는 방법으로 그림 9(b)와 같이 변환, 차분 및 바이어스 등을 활용, 데이터 변환을 통해 분산 정상화를 수행한다. 그 결과, 1초부터 시스템 불안정 전 24초까지 평균 $-1.865\times 10^{-6}$, 표준편차 $0.000456$으로 통계적으로 분석하기 유의한 데이터를 취득할 수 있다.

그리고 AR 모형의 적합성과 계수의 차수를 선정하기 위하여, 그림 10(a)와 같이 자기상관함수 ACF (auto correlation function)와 그림 10(b)와 같이 편자기상관함수 PACF (Partial auto correlation function)를 통한 통계적인 분석을 수행한 뒤, 자기회귀 모형 기반으로 안정도를 추정한다.

그림. 9. $i_{d}$ 전류값에 대한 시계열 정상화 전‧후 비교 (a) 시계열 정상화 전 (b) 시계열 정상화 후.

Fig. 9. Comparison between stationarizing time series data (a) Before stationarizing process, and (b) after stationarizing process.

이때, AR 모형을 적용하기 위한 ACF 및 PACF의 특징을 살펴보면, 기하학적으로 ACF는 지수함수적으로 0으로 tail-off되는 형상을 볼 수 있으며, PACF의 경우 lag 22 지점에서 절단점 (cut-off)을 형성하는 것을 볼 수 있으므로, 본 논문에서 활용하는 $i_{D}$ 전류의 시계열 데이터에 대해 AR 모형을 활용하는 것은 적합하다고 볼 수 있다. 다음 섹션의 안정도 판정을 위한 AR의 차수 또한 Lagging factor 22 샘플로 하여 같게 설정 후 진행하고자 한다.

3.3.1.2 시계열 자기회귀 모형 기반 안정도 추정

기본적으로 자기회귀 모형은 아래 수식과 같이 나타낼 수 있다 ⁽¹⁰⁾, ⁽¹¹⁾.

(8)

$Y_{t}=c+\sum_{i=1}^{p}\varphi_{i}Y_{t-i}+\varepsilon_{t}$

단, $\varphi_{1},\:\cdots ,\:\varphi_{p}$는 AR(p) 시계열 모형의 계수를 의미하고, $c$는 바이어스된 상수, 그리고 $\varepsilon_{t}$는 노이즈를 나타낸다. 여기서 위 계수 $\varphi_{i}$를 추정하는 방법은 Yule-Walker method를 활용하였고, 이는 안정성이 무조건 보장되는 방법이므로 채용하였다. 그리고 적당한 lagging을 둔 AR(22) 모형을 구축하였으므로, 이는 결국 이산시간 시스템의 차분방정식의 형태로 볼 수 있으며, z-영역의 전달함수로 표현할 수 있다. 그리고 연속시간 시스템의 주파수 s-영역으로 나타내기 위하여, Tustin 기법으로 변환하여 나타낼 수 있게 되며, 그 변환 과정은 다음과 같다.

(9)

$z=e^{s T_{s}}=\dfrac{e^{s T_{s}/2}}{e^{-s T_{s}/2}}\approx\dfrac{1+s T_{s}/2}{1-s T_{s}/2}$

$z$는 이산주파수 영역, $s$ 는 연속 주파수 영역의 도메인, $T_{s}$ 시간격을 나타낸다. 이렇게 연속시간 시스템에서 전달함수로

그림. 10. $i_{D}$에 대한 ACF와 PACF 분석 (d축 전류) (a) ACF (자기상관함수)와 (b) PACF (편자기상관함수) 22 샘플의 지연에서 절단점 형성.

Fig. 10. ACF and PACF of stationarized data (a) ACF, and (b) PACF and cut-off point of lagging factor of 22 sample.

표현될 수 있다면, 상태공간 방정식으로 표현하는 것에 무리가 없으므로, 다음과 같은 관계에 의해 미분방정식으로 나타낼 수 있으며, 최종적으로 AR(p) 모형의 지연 요소 p와 같은 차원인 시스템의 고유값을 구할 수 있다.

(10)

$\dfrac{d^{n}y}{dt^{n}}+a_{n-1}\dfrac{d^{n-1}y}{dt^{n-1}}+\cdots +a_{1}\dfrac{dy}{dt}+a_{0}y=b_{0}u$

여기서, $y$와 $u$는 각각 미분방정식의 출력과 입력변수를 나타내며, $a_{0},\:\cdots ,\: a_{n-1}$과 $b_{0}$는 AR(p) 시계열 모형의 미분방정식 중 출력과 입력에 해당되는 계수를 각각 나타낸다.

(11)

$\{\begin{aligned}\dot X =A\bullet X+B\bullet U\\Y = C\bullet X +D\bullet U\end{aligned}$

단, 위 식은 상태공간 모델의 구성 요소로서, $X$는 상태변수 벡터이고, $Y$는 출력벡터, 그리고 $U$는 입력벡터를 나타낸다. 또한, $A$는 상태행렬, $B$는 입력행렬, $C$는 출력행렬, 그리고 $D$는 피드포워드 행렬을 의미한다.

그림 11은 위 과정 및 스킴을 개략적으로 반영한 구성도로 나타낸 것이다. 이때, OPAL-RT의 RT-LAB 및 MATLAB Simulink로 구현함을 고려하여 샘플링 요소, 변수 선택 및 배열형 데이터의 buffer 등을 반영하였다. 또한, 해당 데이터 처리 스킴은 모델 내 Master Subsystem (SM)에 구현하였다. 여기서 AR 뿐만 아니라 추후 실시간 시뮬레이션에 시간 기반 제어 및 감시 요소를 활용하기 위해서 실질적으로 연산을 고려해야 하므로, buffer 블록을 도입한다. 또한, AR(p)도 p차의

그림. 11. PMSM 견인전동기의 실시간 시뮬레이션에서 수치적 안정도 분석을 위한 기본적인 신호처리 프로세스.

Fig. 11. The basic scheme of procedure of the numerical stability analysis for the real-time simulation with PMSM.

지연을 고려해야 하므로, 샘플링 시간에 총 $k$배와 $p$배를 적용한 추정 지연이 생기게 되어 이를 고려해 실시간 시뮬레이션 시간격을 조정해야 한다.

3.4 통계적 기법의 실시간 시뮬레이션 연산에 대한 영향

해당 시뮬레이션은 overrun이 발생하지 않도록 시뮬레이션 시간격 $T_{S}$를 50 $\mu s$로 설정하였다. 일반적으로 스위칭 신호가 많은 철도차량의 경우 더 낮게 설정하는 것이 보통이나, 본 논문에서는 MATLAB에 내장된 SVPWM (Space Vector Voltage Pulse width Modulation) Generator 블록을 스위칭주파수를 10kHz로 낮은 대역에서 사용하여, 주파수와 듀티비만 입력되면 되므로, 삼각파를 생성할 필요가 없다. 또한 위에서 언급하였듯, 사용한 통계적 처리 알고리즘을 활용하기 위해서 $k$배와 $p$배의 지연요소를 고려하기 때문에 길게 설정한 것이다. 이때, 시뮬레이션 성능 (연산 시간 점유율)을 비교할 경우, 기본 구성 모델의 경우 핵심 연산 시간이 평균 4.8$\mu s$ 정도이지만, 취득 데이터 분석 알고리즘을 포함한다면, 평균 5.5 $\mu s$ 이상으로 소폭 상승하는 것으로 나타낸다. 다만, 이것은 같은 시스템 내에 다른 시 간격을 갖고 있는 부분 (sub-rate)의 비중이 작을 경우이고, 해당 모델의 개발 과정 중 내부 신호 제어 및 신호처리를 위하여 다양한 시 간격을 적용하였을 때는 overrun이 다수 관찰된 바 있다. 따라서 복잡한 신호처리 등 분석 기법이 모델 내에 포함되어도 실시간성에 크게 어긋나지 않으나, 모델 내의 신호의 흐름이 원활하도록 시 간격을 잘 맞추어 설정해 주어야 한다.

그림. 12. PMSM 제어시스템에서 안정 영역부터 불안정 영역까지의 동적특성을 나타내는 AR 기반 시스템 안정도 스캐닝 기법. (a) 연속시간시스템의 주파수 평면, (b) 이산시간시스템의 주파수 평면.

Fig. 12. A visualization of AR based system stability scanning from initial stage to unstable stage of PMSM control system (a) Continuous time system frequency plane, (b) Discrete time system frequency plane.