C++를 활용한 블록선도 예제

일반적인 블록선도는 제어 시스템을 거시적으로 관찰하고 분석하는 강력한 도구다. 전달함수들을 직렬 연결, 병렬 연결, 피드백 연결로 조합하여 시스템 응답을 살펴볼 수 있다. 이를 C++ 코드로 구현한다면, 이산화된 알고리즘을 통해 제어시스템의 특성 분석부터 시뮬레이션까지 수행할 수 있다. 특히 오브젝트 지향 패러다임을 활용하여, 제어 블록(전달함수)을 클래스로 구현하고 여러 블록을 중첩하여 복합적인 제어 시스템을 손쉽게 프로그래밍할 수 있다.

블록선도의 이론적 배경

블록선도에서 핵심은 전달함수를 정의하고, 이를 다양한 방식(직렬, 병렬, 피드백)으로 연결하여 전체 전달함수를 유도하는 과정이다. 라플라스 변수 $s$에 대해 시스템의 입력을 $R(s)$, 출력이 $Y(s)$라 하면, 가장 단순한 1차 전달함수는 다음과 같이 표현된다.

$$G(s) = \frac{K}{\tau s + 1}$$

이때 $K$와 $\tau$는 각각 이득(gain)과 시정수(time constant)를 의미한다. 여러 개의 전달함수를 직렬로 연결하면, 전반적인 전달함수는 두 전달함수의 곱으로 주어진다. 가령 $G_1(s)$과 $G_2(s)$가 직렬 연결되어 있다고 할 때, 전체 전달함수는

$$G_{\text{series}}(s) = G_1(s) \times G_2(s)$$

병렬 연결 시에는 두 전달함수가 합으로 나타나며,

$$G_{\text{parallel}}(s) = G_1(s) + G_2(s)$$

피드백 루프가 있는 경우에는 닫힌루프 전달함수 $T(s)$를

$$T(s) = \frac{G(s)}{1 + G(s) H(s)}$$

로 정의한다. 여기서 $H(s)$는 센서나 피드백 요소를 의미한다. 아래 mermaid 블록선도 예시는 단순한 유니티 피드백 구조를 나타낸다.

C++ 프로그래밍 환경 준비

C++로 제어 알고리즘이나 블록선도 구현 예제를 작성하려면 적절한 컴파일러와 프로그래밍 환경이 갖추어져야 한다. 다음과 같이 G++ 또는 clang++ 등을 사용할 수 있다.

$ g++ main.cpp -o main
$ ./main

또한 Makefile을 작성하여 컴파일 과정을 자동화할 수도 있다. 예를 들어, Linux나 macOS 환경에서 아래와 같은 단순 Makefile을 작성해 사용할 수 있다.

all:
    g++ main.cpp -o main

TransferFunction 클래스 설계

전달함수를 클래스로 추상화하면 여러 구성요소와의 연결 연산을 멤버 함수로 정의할 수 있다. 전형적으로 분자와 분모를 다항식(polynomial) 계수로 표현한다. 예를 들어,

$$G(s) = \frac{b_0 + b_1 s + b_2 s^2 + \cdots + b_m s^m}{a_0 + a_1 s + a_2 s^2 + \cdots + a_n s^n}$$

를 나타내기 위해 분자 계수를 \mathbf{b}, 분모 계수를 \mathbf{a}로 두어 관리한다. 다음은 분자와 분모를 std::vector 형으로 저장하는 간단한 예시다.

#include <vector>
#include <iostream>

class TransferFunction {
private:
    std::vector<double> numerator;   // 분자 계수
    std::vector<double> denominator; // 분모 계수

public:
    // 생성자
    TransferFunction(const std::vector<double>& num,
                     const std::vector<double>& den)
        : numerator(num), denominator(den) {}

    // Copy constructor, assignment operator 등은 생략

    // 직렬 연결
    TransferFunction series(const TransferFunction& other) const {
        // 분자 계수끼리의 합성(곱), 분모 계수끼리의 합성(곱)을 수행
        // 실제 구현은 다항식의 곱으로 표현
        std::vector<double> newNum(numerator.size() + other.numerator.size() - 1, 0.0);
        std::vector<double> newDen(denominator.size() + other.denominator.size() - 1, 0.0);

        // 다항식 곱 계산 (분자)
        for (size_t i = 0; i < numerator.size(); i++) {
            for (size_t j = 0; j < other.numerator.size(); j++) {
                newNum[i + j] += numerator[i] * other.numerator[j];
            }
        }

        // 다항식 곱 계산 (분모)
        for (size_t i = 0; i < denominator.size(); i++) {
            for (size_t j = 0; j < other.denominator.size(); j++) {
                newDen[i + j] += denominator[i] * other.denominator[j];
            }
        }

        return TransferFunction(newNum, newDen);
    }

    // 병렬 연결
    // (G_1 + G_2) = (b1/a1 + b2/a2) = (b1*a2 + b2*a1) / (a1*a2)
    TransferFunction parallel(const TransferFunction& other) const {
        std::vector<double> newNum;
        std::vector<double> newDen;

        // a1*a2 계산
        std::vector<double> denMul(denominator.size() + other.denominator.size() - 1, 0.0);
        for (size_t i = 0; i < denominator.size(); i++) {
            for (size_t j = 0; j < other.denominator.size(); j++) {
                denMul[i + j] += denominator[i] * other.denominator[j];
            }
        }
        newDen = denMul;

        // b1*a2 계산
        std::vector<double> b1a2(numerator.size() + other.denominator.size() - 1, 0.0);
        for (size_t i = 0; i < numerator.size(); i++) {
            for (size_t j = 0; j < other.denominator.size(); j++) {
                b1a2[i + j] += numerator[i] * other.denominator[j];
            }
        }

        // b2*a1 계산
        std::vector<double> b2a1(other.numerator.size() + denominator.size() - 1, 0.0);
        for (size_t i = 0; i < other.numerator.size(); i++) {
            for (size_t j = 0; j < denominator.size(); j++) {
                b2a1[i + j] += other.numerator[i] * denominator[j];
            }
        }

        // newNum = b1*a2 + b2*a1
        // 두 다항식 b1a2, b2a1 더하기
        size_t newSize = std::max(b1a2.size(), b2a1.size());
        newNum.resize(newSize, 0.0);
        for (size_t i = 0; i < b1a2.size(); i++) {
            newNum[i] += b1a2[i];
        }
        for (size_t i = 0; i < b2a1.size(); i++) {
            newNum[i] += b2a1[i];
        }

        return TransferFunction(newNum, newDen);
    }

    // 피드백 연결
    // T(s) = G(s) / (1 + G(s)H(s))
    TransferFunction feedback(const TransferFunction& H) const {
        // G(s)*H(s)를 구한 뒤 1 + G(s)*H(s)의 분자/분모를 형성
        TransferFunction GH = this->series(H);

        // 1 + GH = (GH.den + GH.num) / GH.den 형식이 아니라,
        // 1을 GH.den으로 표현하여 다항식을 더해야 한다.
        // T(s) = G / (1 + GH)
        // 분자:  G.num * GH.den
        // 분모:  GH.den + G.num*H.num (교차 항 고려)
        // 정확한 계산을 위해서는 다항식 덧셈/곱셈을 적절히 구현해야 한다.

        // 분자 = G.num * H.den (이미 GH.den)
        // GH.den = G.den * H.den
        // GH.num = G.num * H.num
        // 분모 = GH.den + GH.num
        //        = (G.den * H.den) + (G.num * H.num)

        std::vector<double> GHnum = GH.numerator;
        std::vector<double> GHden = GH.denominator;

        // 분자: this->numerator * H.denominator
        std::vector<double> newNum(this->numerator.size() + H.denominator.size() - 1, 0.0);
        for (size_t i = 0; i < this->numerator.size(); i++) {
            for (size_t j = 0; j < H.denominator.size(); j++) {
                newNum[i + j] += this->numerator[i] * H.denominator[j];
            }
        }

        // 분모: GHden + GHnum
        // GHden = (this->denominator * H.denominator)
        // GHnum = (this->numerator   * H.numerator  )
        size_t maxSize = std::max(GHden.size(), GHnum.size());
        std::vector<double> newDen(maxSize, 0.0);
        for (size_t i = 0; i < GHden.size(); i++) {
            newDen[i] += GHden[i];
        }
        for (size_t i = 0; i < GHnum.size(); i++) {
            newDen[i] += GHnum[i];
        }

        return TransferFunction(newNum, newDen);
    }

    // 전달함수의 분자와 분모를 화면에 출력
    void print() const {
        std::cout << "Numerator: ";
        for (auto c : numerator) {
            std::cout << c << " ";
        }
        std::cout << std::endl;

        std::cout << "Denominator: ";
        for (auto c : denominator) {
            std::cout << c << " ";
        }
        std::cout << std::endl;
    }
};

위 코드에서는 다항식의 곱셈과 덧셈을 엄밀하게 구현하여 직렬, 병렬, 피드백 연결을 모두 표현할 수 있다. 예를 들어, 직렬 연결은 다항식 곱으로, 병렬 연결은 통분을 통한 다항식 덧셈으로 표현된다. 피드백 연결의 경우 $1 + G(s)H(s)$ 부분이 다항식 덧셈으로 처리됨을 볼 수 있다. 실제 제어 알고리즘을 설계할 때는 이 클래스의 멤버 함수를 적절히 조합하여 전체 시스템 전달함수를 단계적으로 구성해 나가면 된다.

TransferFunction 클래스 활용 예시

TransferFunction 클래스를 활용해 간단한 블록선도를 구성해 볼 수 있다. 예를 들어 다음과 같은 1차, 2차 시스템 전달함수들을 고려해 보자.

$$\begin{aligned} G_1(s) &= \frac{1}{s+1} \\ G_2(s) &= \frac{s + 2}{s^2 + 3s + 1} \end{aligned}$$

여기서 $G_1(s)$는 분자(numerator)가 ${1}$, 분모(denominator)가 ${1, 1}$이고, $G_2(s)$는 분자 계수가 ${2, 1}$ (상수항부터 s항 순서), 분모 계수가 ${1, 3, 1}$ (상수항, s항, s^2 항 순서)로 구성될 수 있다.

이 두 전달함수를 C++ 코드로 생성하고, 이를 직렬로 연결한 $G_\text{series}(s) = G_1(s) G_2(s)$, 병렬로 연결한 $G_\text{parallel}(s) = G_1(s) + G_2(s)$, 유니티 피드백으로 감싼 $T(s) = \frac{G_1(s)}{1 + G_1(s)}$ 등등을 계산해 보도록 하자. main.cpp에 아래와 같은 코드를 작성해 볼 수 있다.

#include <iostream>
#include "TransferFunction.hpp" // 위에서 정의한 클래스가 들어있다고 가정

int main() {
    // G1(s) = 1 / (s + 1)
    // 분자: {1}, 분모: {1, 1}
    TransferFunction G1({1.0}, {1.0, 1.0});

    // G2(s) = (s + 2) / (s^2 + 3s + 1)
    // 분자: {2, 1}, 분모: {1, 3, 1}
    TransferFunction G2({2.0, 1.0}, {1.0, 3.0, 1.0});

    // 직렬 연결
    TransferFunction G_series = G1.series(G2);
    std::cout << "=== G1 * G2 (직렬 연결) ===" << std::endl;
    G_series.print();
    std::cout << std::endl;

    // 병렬 연결
    TransferFunction G_parallel = G1.parallel(G2);
    std::cout << "=== G1 + G2 (병렬 연결) ===" << std::endl;
    G_parallel.print();
    std::cout << std::endl;

    // 유니티 피드백 ( H(s) = 1 )
    // T(s) = G1 / (1 + G1 * 1)
    // 여기서 H(s)=1이므로 분자: {1}, 분모: {1}
    TransferFunction H({1.0}, {1.0});
    TransferFunction T = G1.feedback(H);
    std::cout << "=== G1 유니티 피드백 닫힌루프 ===" << std::endl;
    T.print();
    std::cout << std::endl;

    return 0;
}

코드를 컴파일 후 실행하면, $G_1(s)$와 $G_2(s)$를 직렬로 연결한 경우의 분자와 분모, 병렬로 연결한 경우의 분자와 분모, 그리고 유니티 피드백 닫힌루프 전달함수의 분자와 분모를 확인할 수 있다.

실행 결과로 출력되는 수치 값들은 실제로 라플라스 도메인 상에서 $G_1(s) G_2(s)$, $G_1(s) + G_2(s)$, $\frac{G_1(s)}{1+G_1(s)}$에 대응하는 다항식 형태가 된다. 예를 들어, $G_\text{series}(s) = \frac{(1)(s+2)}{(s+1)(s^2+3s+1)} = \frac{s+2}{s^3 + 4s^2 + 4s + 1}$ 같은 식으로 확장되는 과정을 코드에서 내부적으로 처리해 준다.

이 예시에 맞추어 코드를 확장하면, 사용자 정의 전달함수를 입력받아 자동으로 직렬, 병렬, 피드백 연결을 계산하고, 출력 형태 역시 원하는 대로 편집하거나 라플라스 변환의 역변환(부분 분수 분해 등)을 적용하여 시간 응답 형태를 분석하는 방향으로 나아갈 수 있다.

시뮬레이션과 시간 응답 계산 확장

블록선도를 정의하고 전달함수를 연산으로 조합하는 것만으로도 유익하지만, 실제 시스템 동작을 모사하기 위해서는 시간 영역 시뮬레이션도 고려해야 한다. 예컨대

$$G(s) = \frac{Y(s)}{U(s)}$$

형식으로 정의된 시스템에 대해, 입력 $\mathbf{u}(t)$가 주어졌을 때 출력 $\mathbf{y}(t)$를 수치 적분 알고리즘(오일러법, 룽게-쿠타법 등)으로 시뮬레이션할 수 있다. 이 과정에서 $G(s)$를 상태방정식 형태로 변환하거나, 이산화(discretization)하여 $G(z)$로 바꾼 뒤에 구현한다.

TransferFunction 클래스를 그대로 이용해서 시뮬레이션을 수행하려면, 추가적으로 미분 방정식을 풀거나, $Z$-변환 방식으로 변환하여 누적 합 연산을 통해 시뮬레이션을 진행하는 함수를 구현할 수 있다. 이를 위해서는 제어 이론에서 다음과 같은 표준 변환 관계를 고려할 수 있다.

$$y(t) + a_1 \frac{dy(t)}{dt} + \dots + a_n \frac{d^n y(t)}{dt^n} = b_0 u(t) + b_1 \frac{du(t)}{dt} + \dots + b_m \frac{d^m u(t)}{dt^m}$$

이 식을 전진 오일러나 룽게-쿠타 같은 기법으로 수치 해석하면, $t=0$부터 $t=T$까지의 구간에서 $\mathbf{y}(t)$, $\mathbf{u}(t)$ 등을 순차적으로 계산할 수 있다.

C++ 코드로 구현할 때는, TransferFunction 클래스 안에 “시간 응답 계산 함수”를 추가하거나, 외부에서 이 클래스를 받아 처리하는 별도의 Simulate 함수를 작성할 수 있다. 이러한 구조화를 통해 블록 단위 제어 시스템과 시뮬레이션 로직이 분리되므로, 더 편리하게 코드 유지보수가 가능하다.

상태방정식으로의 확장

전달함수 기반의 블록선도 설계는 단입력 단출력(SISO) 시스템에서 직관적이고 간편하게 이루어진다. 하지만 다입력 다출력(MIMO) 구조나 비선형 요소가 포함된 경우, 혹은 시간영역에서의 시뮬레이션이 주가 되는 경우에는 상태방정식(State-Space) 표현이 유리하다. 상태방정식을 사용하면, 고차 시스템이나 복합 시스템에 대해서도 체계적으로 구성된 행렬 연산만으로 시뮬레이션, 설계, 분석이 가능하다.

전형적인 상태방정식은 다음과 같이 표현한다.

$$\dot{\mathbf{x}}(t) = \mathbf{A}\,\mathbf{x}(t) + \mathbf{B}\,\mathbf{u}(t) \\ \mathbf{y}(t) = \mathbf{C}\,\mathbf{x}(t) + \mathbf{D}\,\mathbf{u}(t)$$

여기서

• $\mathbf{x}(t)$는 상태벡터,

• $\mathbf{u}(t)$는 입력벡터,

• $\mathbf{y}(t)$는 출력벡터이며,

• $\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}$는 시스템 행렬이다.

전달함수 $G(s)$가 한 개의 입력과 한 개의 출력을 갖는 SISO 형태라면, 이를 최소실현(minimal realization) 형태의 상태방정식으로 변환하는 여러 기법이 존재한다. 예를 들어 제어공학에서 잘 알려진 방법인 직렬/병렬/캐논형(canonical form)이 그 예다. 차수가 $n$인 전달함수

$$G(s) = \frac{b_0 + b_1 s + \dots + b_m s^m}{a_0 + a_1 s + \dots + a_n s^n}$$

를 고려하면, 이를 다음과 같은 캐논형으로 표현할 수 있다.

가령 제어행렬 $\mathbf{B}$, 관측행렬 $\mathbf{C}$, 직접이득행렬 $\mathbf{D}$ 등을 구성하고, C++ 코드에서 동적 할당으로 $\mathbf{A}$ 행렬을 만든 뒤, $t$를 일정 간격으로 증가시키며 전진 오일러 또는 룽게-쿠타와 같은 수치 적분 방식으로 $\mathbf{x}(t + \Delta t)$를 업데이트하면 된다. 예를 들어 전진 오일러법은

$$\begin{aligned} \mathbf{x}(k+1) &= \mathbf{x}(k) + \Delta t \left(\mathbf{A} \mathbf{x}(k) + \mathbf{B} \mathbf{u}(k)\right) \\ \mathbf{y}(k+1) &= \mathbf{C} \mathbf{x}(k+1) + \mathbf{D} \mathbf{u}(k+1) \end{aligned}$$

형식으로 구현할 수 있다.

간단한 상태방정식 클래스 예시

C++에서 상태방정식을 다루기 위해 간단한 클래스를 만들어, 행렬 곱셈과 벡터 연산을 추상화할 수 있다. 예를 들어 다음과 같은 구조로 짜볼 수 있다.

#include <vector>
#include <iostream>

class StateSpace {
private:
    std::vector<std::vector<double>> A;
    std::vector<std::vector<double>> B;
    std::vector<std::vector<double>> C;
    std::vector<std::vector<double>> D;

    std::vector<double> x; // 현재 상태벡터

public:
    // 생성자
    StateSpace(
        const std::vector<std::vector<double>>& A_in,
        const std::vector<std::vector<double>>& B_in,
        const std::vector<std::vector<double>>& C_in,
        const std::vector<std::vector<double>>& D_in,
        const std::vector<double>& x0
    )
    : A(A_in), B(B_in), C(C_in), D(D_in), x(x0)
    {}

    // 간단한 전진 오일러법에 의한 상태 업데이트
    void updateEuler(const std::vector<double>& u, double dt) {
        // dx = A*x + B*u
        std::vector<double> dx(A.size(), 0.0);

        // dx = A*x
        for (size_t i = 0; i < A.size(); i++) {
            for (size_t j = 0; j < A[0].size(); j++) {
                dx[i] += A[i][j] * x[j];
            }
        }

        // dx += B*u
        for (size_t i = 0; i < B.size(); i++) {
            for (size_t j = 0; j < B[0].size(); j++) {
                dx[i] += B[i][j] * u[j];
            }
        }

        // x = x + dt * dx
        for (size_t i = 0; i < x.size(); i++) {
            x[i] += dt * dx[i];
        }
    }

    // 출력 y = C*x + D*u
    std::vector<double> output(const std::vector<double>& u) const {
        std::vector<double> y(C.size(), 0.0);

        // y = C*x
        for (size_t i = 0; i < C.size(); i++) {
            for (size_t j = 0; j < C[0].size(); j++) {
                y[i] += C[i][j] * x[j];
            }
        }

        // y += D*u
        for (size_t i = 0; i < D.size(); i++) {
            for (size_t j = 0; j < D[0].size(); j++) {
                y[i] += D[i][j] * u[j];
            }
        }

        return y;
    }

    // 현재 상태 반환
    const std::vector<double>& getState() const {
        return x;
    }

    // 상태 초기화
    void reset(const std::vector<double>& x0) {
        x = x0;
    }
};

이와 같은 StateSpace 클래스를 통해 MIMO든 SISO든 임의의 차수를 가지는 선형 시스템을 통합적으로 다룰 수 있다. 또한 TransferFunction으로 표현된 시스템을 특정 캐논형 상태방정식으로 변환하는 함수(예: TF2SS)를 별도로 구현하여, 변환된 계수를 StateSpace 객체로 만드는 기능도 제공할 수 있다.

블록선도와 상태방정식의 결합

블록선도를 상태방정식으로 바꿔서 시뮬레이션하거나, 반대로 상태방정식을 블록선도 형태로 표현하는 일은 복합 시스템 해석에서 자주 쓰인다. 예를 들어 두 개의 SISO 전달함수를 직렬 연결하면, 전체 상태방정식은 두 하위 시스템을 적절히 이어붙인 형태로 조합할 수 있다. 이때 계수 행렬을 블록(block) 단위로 배치하는 기법이 활용된다.

블록선도의 각 요소를 상태방정식으로 변환한 뒤, 전체 시스템을 하나의 큰 상태방정식으로 합쳐 시뮬레이션할 수도 있고, 부분 시스템들 각각을 별도 StateSpace 객체로 선언하여 입력-출력을 연결하는 방법으로 시뮬레이션할 수도 있다. 후자의 방식은 객체지향적이며, 시스템 구성도를 코드 구조와 1:1 매핑하기 쉽다는 장점이 있다.

예컨대 시스템 전체가 아래처럼 직렬 연결된 두 상태공간 $G_1, G_2$라 하면,

$$u \rightarrow [G_1] \rightarrow [G_2] \rightarrow y$$

$G_1$의 출력이 $G_2$의 입력이 되도록, C++ 코드에서 매 시뮬레이션 스텝마다 다음과 같은 순서를 수행할 수 있다.

1. $G_1$에 $u$를 입력해 상태 업데이트 및 출력 $\mathbf{y}_1$ 계산

2. 그 출력 $\mathbf{y}_1$을 $G_2$의 입력으로 하여 상태 업데이트 및 최종 출력 $\mathbf{y}$ 계산

이를 일반화하면 병렬 연결, 피드백 연결, 하위 블록 여러 개를 계층적으로 묶어낸 복합 시스템을 객체 단위로 조립해 나갈 수 있다.

다중입출력 블록선도와 MIMO 상태공간

실제 산업 현장이나 실험실에서 다루는 많은 제어 시스템은 단순 SISO 구조가 아닌, 다중입출력(MIMO) 특성을 갖는다. 이는 하나의 물리계에 여러 개의 입력(제어량)과 여러 개의 출력(측정량)이 존재함을 의미한다. 예컨대 로봇팔의 각 조인트마다 모터 전류를 입력으로 가하고, 각 조인트의 각도와 각속도를 출력으로 측정할 수 있다.

이때 전달함수 기반으로 MIMO 시스템을 표현하면 수많은 전이함수(transfer channel)를 행렬 형태로 정리하게 된다. 즉

$$\mathbf{G}(s) = \begin{bmatrix} G_{11}(s) & \dots & G_{1m}(s) \\ \vdots & \ddots & \vdots \\ G_{n1}(s) & \dots & G_{nm}(s) \end{bmatrix}$$

형식의 행렬전달함수가 구성될 수 있다. 이들을 블록선도로 나타내면 입력 채널별로 각 행 단위의 블록들이 병렬로 존재하고, 행렬의 각 요소 $G_{ij}(s)$는 한 개의 스칼라 전달함수 블록에 해당한다.

하지만 MIMO 시스템의 해석, 시뮬레이션, 그리고 제어기 설계 작업은 상태방정식(Form)에 기반을 두었을 때 훨씬 깔끔해진다. StateSpace 클래스를 사용해 $\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}$ 행렬로 구성해두면, 시스템 차수가 커지더라도 행렬 연산 한 번에 여러 출력 값을 동시에 계산할 수 있기 때문이다. MIMO 상태방정식은 입력벡터와 출력벡터의 차원만 늘어나도 같은 형태를 유지하므로, SISO 시스템과 마찬가지로 전진 오일러, 룽게-쿠타 같은 공통의 수치 적분 기법을 동일하게 적용할 수 있다.

이러한 구조를 C++ 코드로 설계할 때, 행렬 연산 라이브러리(Eigen, Armadillo 등)를 사용하면 구현 효율이 크게 높아진다. StateSpace 클래스에 내부적으로 Eigen::MatrixXd 등을 이용해 $\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}$를 다룬다면, 사용자 정의 연산자를 활용해 빠르고 간결하게 상태방정식을 계산할 수 있다. 예컨대 다음과 같은 패턴으로 코드를 수정할 수 있다.

#include <Eigen/Dense>

class StateSpaceEigen {
private:
    Eigen::MatrixXd A;
    Eigen::MatrixXd B;
    Eigen::MatrixXd C;
    Eigen::MatrixXd D;

    Eigen::VectorXd x; // 상태벡터

public:
    StateSpaceEigen(
        const Eigen::MatrixXd& A_in,
        const Eigen::MatrixXd& B_in,
        const Eigen::MatrixXd& C_in,
        const Eigen::MatrixXd& D_in,
        const Eigen::VectorXd& x0
    )
    : A(A_in), B(B_in), C(C_in), D(D_in), x(x0) {}

    void updateEuler(const Eigen::VectorXd& u, double dt) {
        Eigen::VectorXd dx = A*x + B*u;
        x = x + dt*dx;
    }

    Eigen::VectorXd output(const Eigen::VectorXd& u) const {
        return C*x + D*u;
    }

    const Eigen::VectorXd& getState() const {
        return x;
    }

    void reset(const Eigen::VectorXd& x0) {
        x = x0;
    }
};

이로써 입력과 상태, 출력의 모든 계산이 벡터·행렬 연산자로 표현되어, 시스템 차원과 관계없이 동일한 코드 구조를 유지하게 된다. 로봇, 항공기, 프로세스 제어 시스템 등에서 복수의 모터, 센서가 연결된 고차원 MIMO 시스템을 다룰 때 특히 유용하다.

블록선도의 객체지향적 구성

MIMO 시스템일지라도, 블록선도 관점에서 여전히 여러 개의 상태방정식 객체를 인스턴스로 생성한 뒤, 출력-입력을 서로 연결해 시뮬레이션할 수 있다. 각 객체는 내부적으로 자신만의 $\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}$ 행렬과 상태벡터를 관리하므로, 시스템 전체가 복잡해지더라도 객체 단위로 모듈화가 이루어진다. 예컨대 매 시뮬레이션 스텝에서, 하위 블록부터 순서대로 updateEuler를 호출하면서 출력을 구하고, 상위 블록의 입력으로 전달한 뒤, 최종 출력에 도달하는 식으로 작동한다.

이러한 객체지향적 접근은, 제어 시스템의 구성 요소(플랜트, 컨트롤러, 센서, 외란 모델 등)를 각각의 C++ 클래스 객체로 만든 뒤, 입력-출력 포트를 연결하는 방식으로 자연스럽게 확장된다. 다이어그램을 이용해 가상적인 MIMO 블록선도를 나타내면 다음과 같이 표현할 수 있다.

이는 $2\times2$ 행렬전달함수를 블록선도로 그린 예시다. 실제 코드에서는 $G_{11}, G_{12}, G_{21}, G_{22}$ 각각을 상태방정식 객체 혹은 TransferFunction 객체로 만들어, 유입되는 입력벡터를 적절히 분리하여 updateEuler, output을 호출한 뒤, 최종 결과를 합쳐서 $\mathbf{y}(t)$로 계산할 수 있다.

고차 필터링, 옵서버, 내부 모델 제어

블록선도와 상태공간 모델의 장점을 결합하면, 고차 필터를 설계하여 센서 노이즈를 제거하거나, 루엔버거 옵서버(관측기)를 추가하여 측정 불가능한 상태를 추정하는 일 등이 손쉽게 이루어진다. 이를 코드로 구현할 때, 필터나 옵서버 역시 또 다른 StateSpace 혹은 TransferFunction 블록으로 추가해주면 된다.

예컨대, 센서 출력에 노이즈가 심한 경우, 저역 통과 필터(저주파 필터)를 설계하여 센서 측정값을 후처리할 수 있다. 이때 저역 통과 필터 역시 전달함수 형태

$$F(s) = \frac{\omega_c}{s + \omega_c}$$

로 구현할 수 있으므로, 간단한 1차 시스템 형태의 TransferFunction 혹은 StateSpace로 만들 수 있다. 플랜트 출력이 $y(t)$라면, 실제 측정값 $\hat{y}(t) = F(s),y(t)$로 대체하여 제어기를 안정화할 수 있다.

옵서버(관측기) 설계에 있어서도 비슷한 원리로, 옵서버의 $\mathbf{A}\text{obs}, \mathbf{B}\text{obs}, \mathbf{C}\text{obs}, \mathbf{D}\text{obs}$ 행렬을 상태방정식 객체로 만들어, 측정값과 추정 상태를 순환시키면 된다. 옵서버 이득($\mathbf{L}$)만큼 측정오차를 보정하여 상태를 업데이트하는 형태로, 블록선도 상에서 따로 옵서버 블록을 배치해 입력을 받아 상태를 추정하는 과정을 그대로 코드로 구현하는 셈이다.

이렇듯 C++에서의 블록선도 구현은, 단순한 전달함수 연산부터 고차원 MIMO 상태공간까지 확장할 수 있으며, 그 사이에는 필터링, 옵서버, 내부 모델 제어(IMC) 등 다양한 고급 제어 기법들을 자연스럽게 반영할 수 있다. 모든 것은 블록 단위로 객체화하여 설계할 수 있으므로, 제어 공학적 개념과 코드 구조의 일관성을 유지하기에 용이하다.

고급 제어 기법: LQR, H∞, MPC 등과의 연계

블록선도 설계와 상태방정식을 활용하여 시간영역 시뮬레이션이 가능해지면, 더 나아가 고급 제어 기법을 적용하는 문턱에 이르게 된다. LQR(Linear Quadratic Regulator), H∞ 제어, 모델 예측 제어(MPC) 등 다양한 선형·비선형 기법들은 시스템을 하나의 상태방정식 형태로 정식화한 뒤, 특정 비용함수나 성능지표를 최적화하면서 제어 입력을 구하는 식으로 동작한다.

C++ 기반의 블록선도 환경에서 이러한 고급 제어 알고리즘을 추가적으로 구현할 수 있다. 예컨대, LQR은 다음과 같은 비용함수를 최소화한다.

$$J = \int_{0}^{\infty} \left( \mathbf{x}(t)^{T} \mathbf{Q}\, \mathbf{x}(t) \;+\; \mathbf{u}(t)^{T} \mathbf{R}\, \mathbf{u}(t) \right) dt$$

여기서 $\mathbf{x}(t)$는 상태벡터, $\mathbf{u}(t)$는 제어입력 벡터, $\mathbf{Q}$와 $\mathbf{R}$은 각각 상태와 제어입력에 대한 가중행렬이다. LQR 제어 law는

$$\mathbf{u}(t) = - \mathbf{K}\, \mathbf{x}(t)$$

로 나타나고, 이때 이득행렬 $\mathbf{K}$는 연립 리카티 방정식을 풀어 얻을 수 있다. 선형 시간불변(LTI) 시스템에서 리카티 방정식을 풀면

$$\mathbf{A}^T \mathbf{P} + \mathbf{P}\,\mathbf{A} - \mathbf{P}\,\mathbf{B}\,\mathbf{R}^{-1} \mathbf{B}^T \mathbf{P} + \mathbf{Q} = 0$$

가 성립하고,

$$\mathbf{K} = \mathbf{R}^{-1}\, \mathbf{B}^T \mathbf{P}$$

를 계산함으로써 얻는다. LQR 알고리즘 구현은 대체로 다음과 같은 절차로 진행된다.

1. 시스템 행렬 $\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}$를 정의한다.

2. 가중행렬 $\mathbf{Q}, \mathbf{R}$을 설정한다.

3. 연립 리카티 방정식을 풀어 $\mathbf{P}$ 행렬을 구한다.

4. $\mathbf{K} = \mathbf{R}^{-1}\mathbf{B}^T \mathbf{P}$로 최적이득행렬을 얻는다.

5. 제어입력 $\mathbf{u}(t) = -\mathbf{K},\mathbf{x}(t)$를 블록선도나 상태공간 시뮬레이션 단계에서 적용한다.

위 과정을 C++ 코드로 작성할 때, 리카티 방정식을 푸는 부분은 보통 선형대수 라이브러리(Eigen, LAPACK 등)를 활용한다. 이를 간단한 예시로 살펴보자.

#include <Eigen/Dense>
#include <iostream>

// LQR를 위한 P 행렬 계산 (간단히 Eigen의 콜레스키 분해 기반 알고리즘 사용)
Eigen::MatrixXd solveRiccati(const Eigen::MatrixXd& A,
                             const Eigen::MatrixXd& B,
                             const Eigen::MatrixXd& Q,
                             const Eigen::MatrixXd& R) {
    // 예시용으로, Eigen에서 제공하는 CARE(Continuous-time Algebraic Riccati Equation) 해석 루틴을 직접 사용하거나,
    // 외부 함수를 이용할 수 있다. 여기서는 단순히 '가짜' 구현이라고 가정하자.
    // 실제로는 Numerics. C++ Slicot, control-toolbox 등 다양한 라이브러리를 쓸 수 있음.

    // 가정상, 임의로 P를 Q로 초기화하여 반환 (실제론 제대로 된 알고리즘이 필요)
    Eigen::MatrixXd P = Q;
    return P;
}

int main() {
    // 시스템 차원 예시
    int n = 2; // 상태변수 수
    int m = 1; // 제어입력 수

    // A, B, Q, R 행렬을 정의
    Eigen::MatrixXd A(n, n);
    Eigen::MatrixXd B(n, m);
    Eigen::MatrixXd Q(n, n);
    Eigen::MatrixXd R(m, m);

    A << 0, 1,
         0, 0;
    B << 0,
         1;
    Q << 10,  0,
          0, 10;
    R << 1;

    // Riccati 방정식 해 P 계산
    Eigen::MatrixXd P = solveRiccati(A, B, Q, R);

    // K = R^-1 * B^T * P
    Eigen::MatrixXd K = R.inverse() * B.transpose() * P;
    std::cout << "Calculated LQR Gain K = \n" << K << std::endl;

    // 예: 제어 입력 u = -K*x
    // 시뮬레이션 중에 상태가 x = [x1, x2]^T라면, u = -K*x로 계산

    return 0;
}

이렇게 도출한 LQR 이득행렬 $\mathbf{K}$를 상태방정식 시뮬레이션에 연결하면, 블록선도에서 “컨트롤러 블록”이 LQR 제어로 작동하는 형태가 된다. 그림으로 나타내면 다음과 같이 표현할 수 있다.

H∞ 제어, MPC, 적응 제어(Adaptive Control) 등 다른 제어 기법들도 크게 다르지 않다. 각 기법마다 리카티 방정식 혹은 최적화 문제를 설정하고, 이를 통해 구해진 제어입력을 상태방정식 객체로 전달해주는 식으로 구현할 수 있다. MPC의 경우, 특정 시계열 예측구간에서 비용함수를 매 시점마다 최소화하여 제어입력을 구하는 구조이므로, 보다 복잡한 최적화 솔버가 필요하다. 그래도 큰 틀에서 바라보면 “시스템(플랜트)로부터 상태·출력을 받아, 제어 로직(최적화·리카티 해 등)을 거쳐 다시 플랜트 입력으로 돌려준다”는 블록선도적 사고방식이 적용된다는 점은 동일하다.

실제 구현 시 고려 사항

고급 제어기를 C++로 작성하여 실시간 또는 준실시간 시뮬레이션, 혹은 임베디드 환경으로 이식할 때 고려해야 할 사항들이 있다. 신호 포화와 양방향 제한(saturation, rate limit) 처리를 어디서 어떻게 할지 결정해야 하고, 외란(obstacle)이나 파라미터 불확실성(uncertainty)에 대한 모델링이 정확하지 않을 경우, 제어기 안정성이 어떻게 훼손될 수 있는지도 주의해야 한다.

이런 문제를 처리하기 위해서는 블록선도 상에 포화 블록(입력 출력 값을 일정 범위로 제한), 비선형 블록(마찰, 데드존 등)을 추가하거나, H∞/μ-분석과 같은 강인제어(robust control) 기법을 이용해 모델 불확실성에 대비하는 방식을 고민해야 한다. 역시, 모든 것은 블록 단위로 분리하여 객체로 구현하고, 서로 연결하는 방식으로 확장할 수 있다.

결국, C++를 활용한 블록선도 예제는 간단한 전달함수 연산이나 상태방정식 시뮬레이션을 넘어, 실제 현장의 다양한 요구사항과 고급 제어 알고리즘을 포괄하는 광범위한 엔지니어링 환경으로 확장 가능하다. 제어공학 입문 단계에서는 상대적으로 작은 예제로 시작하되, 점차 복잡한 모델과 고급 기법을 추가해나가며 전체 시스템의 구조와 작동원리를 체득할 수 있다.

실시간 구현과 하드웨어 인터페이스

C++로 블록선도 기반 시뮬레이션을 구현하다 보면, 이를 실제 하드웨어와 연결하여 실시간 제어 또는 준실시간 제어를 수행하고자 하는 요구가 발생한다. 예컨대 DC 모터를 구동하는 마이크로컨트롤러(MCU)나 산업용 PC를 통해, 이론적으로 설계한 블록선도를 그대로 임베디드 시스템 환경에서 실행해 볼 수 있다. 이때 고려해야 할 사항들을 살펴보자.

실시간 성능 확보를 위해서는 블록선도 연산(상태방정식 업데이트, 전달함수 연산, 옵서버, 제어기 등)이 주기적(periodic)으로 동작해야 한다. 일반적인 PC 환경에서는 운영체제 스케줄러의 간섭으로 인해 시간지연이 발생할 수 있다. 따라서 실시간 운영체제(RTOS)나 실시간 커널 패치를 적용한 리눅스(RT-Linux), 혹은 Windows에서의 하드 실시간 지원(정형화된 미들웨어 등)을 사용하여 주기적 태스크를 보장받을 수 있는 환경을 구성하기도 한다.

하드웨어 인터페이스는 ADC(Analog-to-Digital Converter)를 통해 센서 데이터를 읽어오고, DAC(Digital-to-Analog Converter)나 PWM(Pulse Width Modulation) 출력을 통해 액추에이터를 구동한다. 예컨대 DC 모터 제어 시스템에서 제어 입력은 모터 드라이브에 공급되는 PWM Duty Ratio가 될 수 있고, 출력 측정값은 엔코더(Encoder)나 홀 센서(Hall sensor)로부터 읽은 회전 위치·속도일 수 있다. 블록선도 상에서 $u(t)$가 제어기에서 계산된 값이라고 한다면, 실시간 코드에서는 해당 $u(t)$를 PWM 신호로 변환하여 하드웨어 레지스터에 기록하는 과정을 포함해야 한다.

C++ 코드에서 이러한 입출력을 제어하기 위해서는 운영체제별, 하드웨어별로 제공되는 라이브러리나 API를 사용해야 한다. 예를 들어, 임베디드 Linux 환경(예: Raspberry Pi)이라면 /sys/class/gpio, /sys/class/pwm 등의 인터페이스나 WiringPi 같은 라이브러리를 사용할 수 있다. 마이크로컨트롤러(ARM Cortex-M, AVR 등)에서는 벤더가 제공하는 HAL(Hardware Abstraction Layer) API를 호출하여 ADC/PWM 레지스터를 직접 다루는 식으로 구현할 수 있다.

실시간 스레드 구성을 위해서는 C++ 표준 라이브러리의 std::thread, std::mutex, std::condition_variable 등을 사용할 수도 있지만, 정밀한 실시간 주기가 보장되지는 않는다. 더 엄밀한 실시간 환경이 필요하다면, 특정 실시간 프레임워크(RTOS, ROS 2의 실시간 노드, EtherCAT 마스터 등)나 커스텀 스케줄링 코드를 이용하여 타이머 인터럽트 기반 주기를 설정하기도 한다. 내부적으로는 타이머 인터럽트가 발생할 때마다 한 주기(cycle)분량의 블록선도 연산을 수행하는 구조다. 예컨대 일정한 주기 $T_s$에 대해, 다음과 같은 루프가 실행된다고 볼 수 있다.

// (의사코드)
void realtimeControlLoop() {
    auto t_next = std::chrono::steady_clock::now();
    while (running) {
        // 1) 센서 신호 읽기
        // 2) 블록선도 업데이트(상태방정식, 제어기, 옵서버, 필터 등)
        // 3) 제어 입력 출력(PWM 신호 갱신 등)

        t_next += std::chrono::microseconds(static_cast<long>(T_s * 1e6));
        std::this_thread::sleep_until(t_next);
    }
}

PC 환경에서는 std::this_thread::sleep_until가 정확한 주기를 보장하지 못할 수 있으므로, RTOS에서는 타이머 ISR(Interrupt Service Routine) 등으로 좀 더 정교하게 구현한다.

블록선도 자체는 크게 달라지지 않지만, 한 번의 사이클에 센서값을 읽고(Filter, Observer 등의 블록을 통과해 추정·필터링), 제어기를 통해 새로운 입력 $u(t)$를 계산, 그리고 그 $u(t)$를 하드웨어에 적용한다는 프로세스가 반복된다는 점이 실시간 제어의 핵심이다. C++ 클래스 구조로 이를 깔끔하게 정리하려면, 다음과 같은 패턴으로 객체 관계를 구성할 수 있다.

하나의 SystemManager 혹은 ControllerManager 객체가 메인 루프를 돌면서 하위 블록(Plant, Controller, Sensor, Observer 등)을 순차적으로 호출한다. 입력(센서값)과 출력(제어값)을 서로 연결해주고, 그 결과가 다시 다음 루프에서의 초기조건이나 상태가 되어 누적(종속)되는 구조다. 이를 다이어그램으로 표현하면 아래와 비슷해진다.

이 모든 과정에서 가장 중요한 점은 주기 $T_s$가 충분히 짧게 설정되어, 실제 물리계의 동적 변화 속도보다 빠른 레이트로 제어가 수행된다는 것이다. 물리계가 빠르게 변화하는데 제어 주기가 너무 크면, 이산적인 제어와 실제 연속 시스템이 부정합을 일으켜 오차가 커진다. 또한 $T_s$를 지나치게 줄이면, 하드웨어 인터페이스나 CPU 연산이 과부하되거나 잡음(sensing noise)에 민감해질 수 있다. 따라서 제어 이론에서 다루는 표본화 주기 선택 기준, 나이퀴스트 샘플링 레이트, 유한차분 시뮬레이션 안정성 등을 종합적으로 고려하여 적절한 $T_s$를 결정한다.

하드웨어 인터페이스와 연동해 블록선도를 구동하게 되면, 순수 소프트웨어 시뮬레이션을 넘어 실제 장비와 상호작용하는 테스트베드가 마련된다. 즉, C++ 코드 상에서 설계한 블록선도를 수정할 때마다 물리계의 응답을 직관적으로 확인할 수 있다. 이는 실험실 수준에서 프로토타이핑을 하거나, 산업 현장에서 소규모 자동화를 구현할 때 매우 빠르고 유연한 접근법이다.

이처럼 C++로 작성된 블록선도 예제는, 충분한 객체지향 설계와 실시간 하드웨어 인터페이스에 대한 배려가 있다면, 비교적 간단한 시뮬레이션 코드에서 출발하여 곧바로 임베디드 제어 시스템으로 이행할 수 있다는 장점이 있다. 특정 상용 제어 소프트웨어(예: Matlab/Simulink) 없이도 순수 C++ 환경에서 제어 알고리즘을 설계, 시뮬레이션, 실시간 구현까지 일관되게 진행 가능하다.

--- 대신 마무리를 위한 요약

지금까지 전달함수 기반의 블록선도에서 출발해, 상태방정식으로 확장하고, 고급 제어 기법(LQR 등)과 시뮬레이션, 그리고 실시간 하드웨어 구현까지 폭넓게 살펴보았다. 핵심 요지를 정리해 보면 다음과 같다.

• 전달함수 블록: 간단한 SISO 시스템을 대상으로, 직렬·병렬·피드백 연결을 다룬다. 다항식의 곱과 통분 등을 통해 전달함수를 손쉽게 합성할 수 있다.

• C++ 클래스 설계: TransferFunction 클래스를 작성해 다항식 연산을 처리하고, 직렬·병렬·피드백 연결 연산을 메서드로 구현한다. 이를 통해 블록선도 연산을 코드로 추상화할 수 있다.

• 상태방정식: 고차 시스템, MIMO 시스템, 고급 제어 기법을 다룰 때는 상태방정식이 더 직관적이고 효율적이다. C++ 코드에서 행렬 연산 라이브러리(Eigen 등)를 활용해 StateSpace 클래스를 구성하고, 전진 오일러·룽게-쿠타로 시간응답을 시뮬레이션할 수 있다.

• 블록선도 객체지향화: 여러 개의 StateSpace 객체(또는 TransferFunction 객체)를 연결해, 블록선도와 유사한 방식으로 구성할 수 있다. 이는 복잡한 MIMO 제어 시스템, 필터, 옵서버 등을 블록 단위로 유지보수하기 쉽다.

• 고급 제어 기법: LQR, H∞, MPC 등 다양한 기법에서 모두 상태방정식을 이용해 제어법칙(이득행렬 등)을 구하고, 이를 시뮬레이션에 적용한다. C++에서 연립 리카티 방정식, 실시간 최적화 솔버 등을 이용해 구현한다.

• 실시간 하드웨어 인터페이스: PC 시뮬레이션을 넘어, 임베디드 환경(마이크로컨트롤러, RTOS, 또는 RT-Linux 등)에서 ADC/PWM 같은 실제 하드웨어 입출력과 연계하여 순수 C++ 코드로 실시간 제어 시스템을 구현할 수 있다. 이 과정에서 주기적 태스크 스케줄링, 샘플링 주기 설계, 신호 포화 처리 등 현실적인 이슈를 해결해야 한다.

이처럼 C++를 활용한 블록선도 예제는 단순 이론 시뮬레이션부터 고급 제어, 실제 하드웨어 연결까지 아우르는 강력한 도구로 발전할 수 있다. 각각의 단계에서는 제어공학 이론(전달함수, 상태방정식, 신호처리)을 충실히 반영하면서, C++ 언어가 제공하는 객체지향 설계, 고성능 연산, 하드웨어 접근성을 적극 활용한다.