문제 정의

전력계통 환경경제 스케줄링(EED)은 발전 비용 최소화와 환경 오염물질 배출량 제한을 동시에 만족해야 하는 대표적인 다목표 최적화 문제다. 본 구현에서는 6대 발전기를 대상으로 연료비와 배출량을 동시에 최적화한다.

알고리 핵심

기본 DE(Differential Evolution)에 비지배 정렬과 혼잡도 거리를 결합한 다목표 최적화 전략을 적용했다. 변이율 F와 교차율 CR을 동적으로 조정하는 적응형 파라미터 메커니즘을 도입해 수렴 속도와 다양성 균형을 맞춘다.

데이터 구조

발전기 특성 데이터는 다음과 같이 구성된다:

• 연료비 계수: 2차 함수 계수 a, b, c 및 밸브점 효과 항
• 배출량 계수: 지수 함수 및 2차 함수 계수
• 출력 한계: 최소/최대 발전 출력
• 손실 렬: B 계수를 통한 전력 손실 계산

핵심 코드 구현

모델 초기화 및 목표 함수

import numpy as np
import matplotlib.pyplot as plt


def parse_system_data(filepath):
    """발전기 데이터 파일 파싱"""
    lines = open(filepath, 'r').readlines()
    gen_count = int(lines[0].split()[-1])
    
    fuel_coeffs = np.zeros((gen_count, 5))
    emission_coeffs = np.zeros((gen_count, 5))
    power_bounds = np.zeros((gen_count, 2))
    loss_matrix = np.zeros((gen_count, gen_count))
    
    # 연료비 계수 및 출력 한계
    for idx in range(gen_count):
        tokens = lines[2 + idx].split()[1:]
        fuel_coeffs[idx, :5] = [float(x) for x in tokens[:5]]
        power_bounds[idx] = [float(tokens[-2]), float(tokens[-1])]
    
    # 배출량 계수
    for idx in range(gen_count):
        tokens = lines[3 + gen_count + idx].split()[1:]
        emission_coeffs[idx, :len(tokens)] = [float(x) for x in tokens]
    
    # 손실 렬 B
    for idx in range(gen_count):
        loss_matrix[idx] = np.array(
            list(map(float, lines[4 + 2*gen_count + idx].split()))
        )
    
    b_vector = np.array(list(map(float, lines[5 + 3*gen_count].split())))
    b_scalar = float(lines[7 + 3*gen_count])
    
    return gen_count, fuel_coeffs, emission_coeffs, power_bounds, \
           loss_matrix, b_vector, b_scalar


def compute_fuel_cost(unit_id, output, coeffs, limits=None):
    """연료비 계산 (밸브점 효과 포함)"""
    base = coeffs[unit_id, 2] * output**2 + coeffs[unit_id, 1] * output + coeffs[unit_id, 0]
    if limits is not None:
        valve_effect = abs(coeffs[unit_id, 3] * 
                          np.sin(coeffs[unit_id, 4] * (limits[unit_id, 0] - output)))
        return base + valve_effect
    return base


def compute_emission(unit_id, output, coeffs, use_expon=True):
    """배출량 계산"""
    if coeffs[0, 3] != 0 and use_expon:
        return (coeffs[unit_id, 0] + 
                (coeffs[unit_id, 1] + coeffs[unit_id, 2]*output) * output +
                coeffs[unit_id, 3] * np.exp(coeffs[unit_id, 4] * output))
    return coeffs[unit_id, 2] * output**2 + coeffs[unit_id, 1] * output + coeffs[unit_id, 0]


class PowerSystemModel:
    def __init__(self, filepath):
        self.n_gen, self.FC, self.EC, self.P_lim, \
        self.B_mat, self.B_vec, self.B_00 = parse_system_data(filepath)
    
    def check_constraints(self, outputs, demand):
        """등식/부등식 제약 조건 검증"""
        total_gen = np.sum(outputs)
        p_loss = np.dot(np.dot(outputs, self.B_mat), outputs)
        p_loss += np.dot(self.B_vec, outputs) + self.B_00
        
        # 전력 균형: 총발전량 = 수요 + 손실
        balance_violation = abs(total_gen - demand - p_loss)
        
        # 출력 한위 위반
        lower_viol = np.sum(np.maximum(0, self.P_lim[:, 0] - outputs))
        upper_viol = np.sum(np.maximum(0, outputs - self.P_lim[:, 1]))
        
        return balance_violation + lower_viol + upper_viol

향상된 다목표 DE 알고리즘

class EnhancedMODE:
    def __init__(self, model, n_iter=1000, pop_size=200, archive_size=100,
                 F_init=0.6, CR_init=0.9):
        self.model = model
        self.n_iter = n_iter
        self.n_pop = pop_size
        self.arc_size = archive_size
        self.F = F_init
        self.CR = CR_init
        
        self.population = None
        self.objectives = None  # [fuel_cost, emission]
        self.archive = []
        self.best_fuel_trace = []
        self.best_emission_trace = []
    
    def initialize_population(self):
        """출력 한위 내에서 균등 분포로 초기화"""
        n = self.model.n_gen
        lower = self.model.P_lim[:, 0]
        upper = self.model.P_lim[:, 1]
        
        self.population = np.random.uniform(
            lower, upper, (self.n_pop, n)
        )
        self.evaluate_all()
    
    def evaluate_all(self):
        """모든 개체의 목표 함수값 계산"""
        self.objectives = np.zeros((self.n_pop, 2))
        for i, individual in enumerate(self.population):
            self.objectives[i] = self._calculate_objectives(individual)
    
    def _calculate_objectives(self, outputs):
        """단일 개체의 연료비와 배출량 계산"""
        fuel = sum(compute_fuel_cost(j, outputs[j], self.model.FC, self.model.P_lim)
                   for j in range(self.model.n_gen))
        emis = sum(compute_emission(j, outputs[j], self.model.EC)
                   for j in range(self.model.n_gen))
        return fuel, emis
    
    def adaptive_mutation(self, target_idx):
        """적응형 차분 변이: F를 동적으로 조정"""
        candidates = [i for i in range(self.n_pop) if i != target_idx]
        a, b, c = np.random.choice(candidates, 3, replace=False)
        
        # 후기 반복에서는 F 감소 (세밀 탐색)
        adaptive_F = self.F * (1 - 0.5 * np.random.rand())
        
        mutant = self.population[a] + adaptive_F * (self.population[b] - self.population[c])
        return np.clip(mutant, self.model.P_lim[:, 0], self.model.P_lim[:, 1])
    
    def crossover(self, target, mutant):
        """교차 연산"""
        mask = np.random.rand(self.model.n_gen) < self.CR
        # 최소 하나는 변이 유전자 포함
        if not np.any(mask):
            mask[np.random.randint(self.model.n_gen)] = True
        return np.where(mask, mutant, target)
    
    def non_dominated_sort(self, objs):
        """비지배 정렬 (NSGA-II 방식)"""
        n = len(objs)
        domination_count = np.zeros(n)
        dominated_sets = [[] for _ in range(n)]
        fronts = [[]]
        
        for i in range(n):
            for j in range(i+1, n):
                if np.all(objs[i] <= objs[j]) and np.any(objs[i] < objs[j]):
                    dominated_sets[i].append(j)
                    domination_count[j] += 1
                elif np.all(objs[j] <= objs[i]) and np.any(objs[j] < objs[i]):
                    dominated_sets[j].append(i)
                    domination_count[i] += 1
            
            if domination_count[i] == 0:
                fronts[0].append(i)
        
        i = 0
        while fronts[i]:
            next_front = []
            for p in fronts[i]:
                for q in dominated_sets[p]:
                    domination_count[q] -= 1
                    if domination_count[q] == 0:
                        next_front.append(q)
            i += 1
            fronts.append(next_front)
        
        return fronts[:-1]  # 마지막 빈 리스트 제외
    
    def crowding_distance(self, objs, front):
        """혼잡도 거리 계산"""
        if len(front) <= 2:
            return {idx: float('inf') for idx in front}
        
        distances = {idx: 0 for idx in front}
        for m in range(2):  # 두 목표 함수에 대해
            sorted_front = sorted(front, key=lambda x: objs[x][m])
            distances[sorted_front[0]] = distances[sorted_front[-1]] = float('inf')
            
            f_min, f_max = objs[sorted_front[0]][m], objs[sorted_front[-1]][m]
            if f_max > f_min:
                for i in range(1, len(sorted_front)-1):
                    distances[sorted_front[i]] += \
                        (objs[sorted_front[i+1]][m] - objs[sorted_front[i-1]][m]) / (f_max - f_min)
        
        return distances
    
    def environmental_selection(self, combined_pop, combined_objs):
        """환경 선택: 비지배 정렬 + 혼잡도 거리"""
        fronts = self.non_dominated_sort(combined_objs)
        selected = []
        
        for front in fronts:
            if len(selected) + len(front) <= self.n_pop:
                selected.extend(front)
            else:
                # 혼잡도 거리로 순위 결정
                distances = self.crowding_distance(combined_objs, front)
                sorted_by_distance = sorted(front, key=lambda x: distances[x], reverse=True)
                selected.extend(sorted_by_distance[:self.n_pop - len(selected)])
                break
        
        self.population = combined_pop[selected]
        self.objectives = combined_objs[selected]
    
    def solve(self, demand):
        """메인 최적화 루프"""
        self.initialize_population()
        
        for generation in range(self.n_iter):
            offspring = []
            off_objs = []
            
            for i in range(self.n_pop):
                # 변이 및 교차
                mutant = self.adaptive_mutation(i)
                trial = self.crossover(self.population[i], mutant)
                
                # 제약 조건 처리: 균등화 보정
                trial = self._repair_solution(trial, demand)
                
                # 목표 함수 평가
                obj = self._calculate_objectives(trial)
                offspring.append(trial)
                off_objs.append(obj)
            
            # 비지배 비교 및 환경 선택
            combined_pop = np.vstack([self.population, np.array(offspring)])
            combined_objs = np.vstack([self.objectives, np.array(off_objs)])
            self.environmental_selection(combined_pop, combined_objs)
            
            # 최적값 기록
            self.best_fuel_trace.append(np.min(self.objectives[:, 0]))
            self.best_emission_trace.append(np.min(self.objectives[:, 1]))
        
        return self.population, self.objectives
    
    def _repair_solution(self, outputs, demand):
        """제약 위반 시 출력 균등화 보정"""
        # 간단한 보정: 총합이 수요에 맞도록 비율 조정
        total = np.sum(outputs)
        if total > 0:
            scale = demand / total
            outputs = outputs * scale
        
        # 한위 클리핑
        outputs = np.clip(outputs, self.model.P_lim[:, 0], self.model.P_lim[:, 1])
        return outputs

실행 및 시각화

if __name__ == '__main__':
    # 시스템 초기화
    system = PowerSystemModel('dataset.txt')
    load_demand = 2.834  # p.u.
    
    # 알고리즘 실행
    optimizer = EnhancedMODE(
        model=system,
        n_iter=1000,
        pop_size=200,
        F_init=0.6,
        CR_init=0.9
    )
    
    solutions, objectives = optimizer.solve(load_demand)
    
    # 수렴 곡선 시각화
    fig, ax1 = plt.subplots(figsize=(10, 6))
    
    color1 = '#E74C3C'
    ax1.set_xlabel('Iteration', fontsize=12)
    ax1.set_ylabel('Minimum Fuel Cost ($/h)', color=color1, fontsize=12)
    line1 = ax1.plot(optimizer.best_fuel_trace, color=color1, linewidth=1.5, label='Fuel Cost')
    ax1.tick_params(axis='y', labelcolor=color1)
    
    ax2 = ax1.twinx()
    color2 = '#3498DB'
    ax2.set_ylabel('Minimum Emission (ton/h)', color=color2, fontsize=12)
    line2 = ax2.plot(optimizer.best_emission_trace, color=color2, linewidth=1.5, label='Emission')
    ax2.tick_params(axis='y', labelcolor=color2)
    
    plt.title('Convergence of Enhanced MODE for EED Problem', fontsize=14)
    fig.tight_layout()
    plt.savefig('convergence_curve.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    # Pareto 전면 출력
    plt.figure(figsize=(8, 6))
    plt.scatter(objectives[:, 0], objectives[:, 1], c='green', alpha=0.6, s=50)
    plt.xlabel('Fuel Cost ($/h)', fontsize=12)
    plt.ylabel('Emission (ton/h)', fontsize=12)
    plt.title('Obtained Pareto Front', fontsize=14)
    plt.grid(True, alpha=0.3)
    plt.savefig('pareto_front.png', dpi=300, bbox_inches='tight')
    plt.show()

확장 적용

본 알고리즘 프레임워크는 전력 스케줄링 외에도 다음 분야로 확장 가능하다:

• 재생에너지 통합 스케줄링: 풍력/태양광 변동성을 고려한 확률적 최적화
• 다지역 경제성장 모델: 지역별 에너지 자립도와 경제성 동시 고려
• 수소 인프라 설계: 생산-저장-운송 비용과 탄소 배출량 최소화

핵심은 목표 함수와 제약 조건의 구조만 교체하면 동일한 다목표 진화 프레임워크를 재사용할 수 있다는 점이다.