군집 지능에서 찾은 최적화의 해답

자연계에서 무리를 지어 이동하는 새 떼는 특별한 지휘자 없이도 효율적으로 먹이를 찾아냅니다. 개체들은 자신의 경험과 주변 동료들의 정보를 바탕으로 이동 방향을 조정하며, 이러한 단순한 상호작용이 전체 집단의 지능적인 행동을 이끌어냅니다. 1995년 Kennedy와 Eberhart는 이 생물학적 메커니즘을 수학적 모델로 형식화하여 입자 군집 최적화(Particle Swarm Optimization, PSO) 알고리즘을 제안했습니다.

초기에는 사회 심리학적 행동 모델링에서 출발했으나, PSO는 함수 최적화, 파라미터 추정, 신경망 학습 등 다양한 분야에서 활용되는 강력한 메타휴리스틱 기법으로 자리 잡았습니다. 본 글에서는 알고리즘의 작동 원리를 분석하고 MATLAB 환경에서의 구현 방법을 단계적으로 살펴봅니다.

알고리즘의 수학적 구조

PSO는 검색 공간 내에 분포된 다수의 후보해를 입자로 모델링합니다. 각 입자는 두 가지 메모리를 유지합니다. 첫째는 자신이 방문했던 위치 중 가장 우수한 해인 개체 최적해(personal best)이고, 둘째는 군집 전체가 발견한 최적해인 전역 최적해(global best)입니다.

매 반복마다 각 입자는 자신의 현재 속도, 개체 최적해로 향하는 인지적 성분, 그리고 전역 최적해로 향하는 사회적 성분을 가중 합산하여 새로운 속도를 계산합니다. 이후 계산된 속도를 바탕으로 위치를 갱신합니다. 속도와 위치 갱신 공식은 다음과 같습니다.

$$ \mathbf{v}_i(t+1) = \omega \mathbf{v}_i(t) + c_1 r_1 (\mathbf{p}_i - \mathbf{x}_i(t)) + c_2 r_2 (\mathbf{g} - \mathbf{x}_i(t)) $$

$$ \mathbf{x}_i(t+1) = \mathbf{x}_i(t) + \mathbf{v}_i(t+1) $$

여기서 $\mathbf{x}_i$와 $\mathbf{v}_i$는 각각 $i$번째 입자의 위치와 속도 벡터, $\mathbf{p}_i$는 개체 최적해, $\mathbf{g}$는 전역 최적해를 의미합니다. $\omega$는 관성 가중치로 이전 속도의 영향력을 조절하며, $c_1$과 $c_2$는 각각 인지적·사회적 학습 인자입니다. $r_1$과 $r_2$는 [0,1] 구간의 난수로 탐색 과정의 확률적 변동을 부여합니다.

MATLAB 구현: 검색 집단 초기화

구현의 첫 단계는 입자 군집을 생성하는 것입니다. 입자의 위치는 검색 공간 내에서 무작위로 배치하고, 속도는 검색 범위에 비례하여 제한된 값으로 초기화합니다. 초기 속도가 너무 크면 입자가 경계를 벗어나며 불안정해지고, 반대로 너무 작으면 수렴 속도가 지나치게 느려집니다.

% 검색 집단 초기화
numParticles = 40;        % 입자 개체 수
numVars = 10;             % 문제 차원
lowerBound = -5 * ones(1, numVars);
upperBound = 5 * ones(1, numVars);

% 입자 위치 생성 (균등 분포)
agentPos = lowerBound + rand(numParticles, numVars) .* (upperBound - lowerBound);

% 속도 한계 설정 (탐색 범위의 15%)
maxVelocity = 0.15 * (upperBound - lowerBound);
agentVel = -maxVelocity + 2 * maxVelocity .* rand(numParticles, numVars);

% 개체 최적해 초기값 설정
indivBestPos = agentPos;
indivBestVal = zeros(numParticles, 1);
for idx = 1:numParticles
    indivBestVal(idx) = evaluateObjective(agentPos(idx, :));
end

% 전역 최적해 탐색
[globalBestVal, bestIdx] = min(indivBestVal);
globalBestPos = agentPos(bestIdx, :);

반복 최적화 메인 루프

초기화 이후에는 반복적으로 적합도 평가, 최적해 갱신, 속도·위치 업데이트를 수행합니다. 이 과정은 종료 조건이 충족될 때까지 지속됩니다. 코드의 가독성과 유지보수성을 위해 각 기능을 별도의 함수로 분리하는 것이 좋습니다.

% PSO 메인 반복
maxCycles = 500;
tolerance = 1e-6;
stagnationLimit = 30;
noImprovement = 0;
prevBest = inf;

for cycle = 1:maxCycles
    % 적합도 산출
    currentVal = evaluateObjective(agentPos);
    
    % 개체 최적해 갱신
    improved = currentVal < indivBestVal;
    indivBestVal(improved) = currentVal(improved);
    indivBestPos(improved, :) = agentPos(improved, :);
    
    % 전역 최적해 갱신
    [cycleMin, minIdx] = min(indivBestVal);
    if cycleMin < globalBestVal
        globalBestVal = cycleMin;
        globalBestPos = indivBestPos(minIdx, :);
    end
    
    % 속도 및 위치 갱신
    agentVel = updateSwarmVelocity(agentPos, agentVel, ...
        indivBestPos, globalBestPos, cycle, maxCycles);
    agentPos = updateSwarmPosition(agentPos, agentVel, lowerBound, upperBound);
    
    % 정체 감지 (수렴 판단)
    if abs(prevBest - globalBestVal) < tolerance
        noImprovement = noImprovement + 1;
    else
        noImprovement = 0;
    end
    prevBest = globalBestVal;
    
    if noImprovement >= stagnationLimit
        fprintf('Early stop at cycle %d (stagnation)\n', cycle);
        break;
    end
end

종료 조건은 최대 반복 횟수뿐만 아니라 연속 정체 감지, 목표 정밀도 달성 등 다양하게 구성할 수 있습니다. 위 코드에서는 일정 횟수 동안 최적해 개선이 없으면 조기 종료하는 로직을 포함했습니다.

속도·위치 갱신 함수의 분리

속도 갱신 공식을 별도 함수로 구현하면 가중치 스케줄링이나 학습 인자 조정 같은 실험을 쉽게 수행할 수 있습니다. 관성 가중치는 반복 진행에 따라 선형적으로 감소시키는 방식을 흔히 사용합니다.

function newVel = updateSwarmVelocity(pos, vel, pBest, gBest, iter, maxIter)
    % 파라미터 설정
    wMax = 0.9;   % 초기 관성 가중치
    wMin = 0.4;    % 종료 시 관성 가중치
    c1 = 1.5;      % 인지적 학습 인자
    c2 = 2.5;      % 사회적 학습 인자
    
    % 선형 감소 관성 가중치
    inertia = wMax - ((wMax - wMin) * iter / maxIter);
    
    % 난수 생성
    r1 = rand(size(pos));
    r2 = rand(size(pos));
    
    % 속도 갱신
    cognitive = c1 * r1 .* (pBest - pos);
    social = c2 * r2 .* (gBest - pos);
    newVel = inertia * vel + cognitive + social;
end

function newPos = updateSwarmPosition(pos, vel, lb, ub)
    newPos = pos + vel;
    
    % 경계 위반 처리 (반사법)
    overLimit = newPos > ub;
    newPos(overLimit) = 2 * ub(overLimit) - newPos(overLimit);
    
    underLimit = newPos < lb;
    newPos(underLimit) = 2 * lb(underLimit) - newPos(underLimit);
    
    % 경계 위반 시 속도 방향 반전 (호출 측에서 처리 필요)
end

경계 처리 전략의 선택

갱신된 위치가 탐색 영역을 벗어날 때 어떻게 처리하느냐에 따라 알고리즘 성능이 크게 달라집니다. 단순히 경계값으로 잘라내는 절단법은 구현이 쉽지만 입자가 벽에 눌려 탐색 다양성이 떨어지는 문제가 있습니다. 반면 경계에 도달하면 튕겨나가는 반사법은 운동 에너지를 보존하면서 탐색을 지속할 수 있게 합니다.

파라미터 튜닝의 기준

PSO 성능은 설정값에 크게 의존합니다. 특히 관성 가중치와 학습 인자의 조합이 탐색-활용 균형을 좌우합니다.

관성 가중치 스케줄링

고정값을 사용하는 것보다 초기에는 크게 설정하여 폭넓은 탐색을 수행하고, 후반부에는 감소시켜 정밀 개선에 집중하는 것이 효과적입니다. 선형 감소 외에도 지수 감소나 비선형 스케줄을 적용할 수 있으며, 다봉 함수에서는 지수 감소가 국소해 함정에서 빠져나오는 데 유리합니다.

% 선형 감소
inertia = 0.9 - (0.5 * cycle / maxCycles);

% 지수 감소
inertia = 0.4 + 0.5 * exp(-3 * cycle / maxCycles);

학습 인자의 비대칭 설정

전통적으로 $c_1 = c_2 = 2.0$을 사용하지만, 사회적 성분을 더 강조하는 $c_1 < c_2$ 설정이 복잡한 지형에서 더 나은 결과를 내는 경향이 있습니다. 다만 지나게 편중하면 조기 수렴 위험이 있으므로 적절한 균형이 필요합니다.

군집 규모 산정

입자 수가 많을수록 탐색 범위는 넓어지지만 계산 비용도 선형적으로 증가합니다. 10차원 이하의 문제에서는 20~30개, 10~30차원에서는 30~50개 정도가 적절하며, 그 이상의 고차원에서는 차원당 1.5~2.5개 비율을 권장합니다.

제약 조건을 포함한 최적화

현실의 문제에는 변수 부호 제한, 합산 한계, 물리적 법칙 등 다양한 제약이 존재합니다. 이를 다루기 위해 벌점 함수법을 활용할 수 있습니다. 제약을 위반하는 정도에 비례하여 적합도에 가중치를 부여해 입자를 자연스럽게 실행 가능 영역으로 유도하는 방식입니다.

function penalizedVal = constrainedObjective(x)
    % 기본 목적 함수
    baseCost = sum(x.^2);
    
    % 제약 위반 패널티 (예: 모든 원소 양수 조건)
    violation = sum(abs(x(x < 0)) .* 1000);
    
    penalizedVal = baseCost + violation;
end

성능 모니터링과 가속화

수렴 과정을 추적하기 위해 매 반복 최적해를 기록하고, 일정 주기마다 입자 분포를 시각화하면 알고리즘 거동을 직관적으로 파악할 수 있습니다. 또한 MATLAB에서는 for 루프보다 행렬 연산 기반의 벡터화가 훨씬 빠르므로, 목적 함수가 벡터 입력을 지원하도록 설계하면 실행 속도를 크게 향상할 수 있습니다.

% 기록용 배열
optimalHistory = zeros(maxCycles, 1);
optimalHistory(cycle) = globalBestVal;

% 분포 시각화 (2차원 투영)
if mod(cycle, 50) == 0
    scatter(agentPos(:,1), agentPos(:,2), '.');
    hold on;
    plot(globalBestPos(1), globalBestPos(2), 'rp', 'MarkerSize', 15);
    title(sprintf('Cycle %d: Opt = %.2e', cycle, globalBestVal));
    drawnow;
end

% 벡터화된 적합도 계산 (Sphere 함수 예시)
currentVal = sum(agentPos.^2, 2);