编写一个程序求解以下多目标优化问题：将n个工件分配给f个工厂，每个工件安排好工厂后只能在该工厂加工，直至加工完成；每个工厂有相同的i个加工阶段，工件在不同的阶段加工时间不同；每个加工阶段有m个并行机器，机器在开始工作后直到该机器上不再安排加工工件才能关闭（即机器在开启后有空闲和加工两种状态，直至机器关闭），机器在各阶段的加工速度不同，机器在各阶段加工时间等于该阶段工件加工时间除以该阶段机器加工速度，每台机器在工作时和空闲时有不同的能耗，总能耗为机器加工时间与加工能耗的乘积加上机器空闲时间与空闲能；总目标是最小化完工时间和最小化总能耗。

这是一个典型的多目标优化问题，可以使用多目标进化算法（例如NSGA-II、MOEA/D等）来求解。以下是一个简单的Python程序，使用DEAP库实现NSGA-II算法来求解该问题：

首先，定义问题的输入参数和目标函数：

import random
from deap import base, creator, tools

# 输入参数
n = 10  # 工件数量
f = 3  # 工厂数量
i = 2  # 加工阶段数量
m = 2  # 并行机器数量

# 加工时间和速度矩阵
proc_time = [[random.uniform(1, 10) for _ in range(i)] for _ in range(n)]
speed = [[random.uniform(1, 10) for _ in range(i)] for _ in range(m)]

# 能耗矩阵
work_energy = [[random.uniform(1, 5) for _ in range(i)] for _ in range(m)]
idle_energy = [[random.uniform(1, 5) for _ in range(i)] for _ in range(m)]

# 定义目标函数
def makespan(individual):
    # 计算加工时间矩阵
    job_machine = [[0 for _ in range(i)] for _ in range(n)]
    for j, f in enumerate(individual):
        for k in range(i):
            machine = min(range(m), key=lambda x: job_machine[j][k][x])
            job_machine[j][k][machine] += proc_time[j][k] / speed[machine][k]

    # 计算完工时间
    return max(max(job_machine[j][i-1]) for j in range(n)),

def energy(individual):
    # 计算加工时间矩阵和能耗矩阵
    job_machine = [[0 for _ in range(i)] for _ in range(n)]
    machine_energy = [[0 for _ in range(i)] for _ in range(m)]
    for j, f in enumerate(individual):
        for k in range(i):
            machine = min(range(m), key=lambda x: job_machine[j][k][x])
            time = proc_time[j][k] / speed[machine][k]
            job_machine[j][k][machine] += time
            machine_energy[machine][k] += time * work_energy[machine][k]

    # 计算空闲时间和能耗
    idle_time = [[max(job_machine[j][k]) - sum(job_machine[j][k][i] for i in range(m) if i != machine)
                  for k in range(i)] for j in range(n)]
    idle_energy = [[time * idle_energy[machine][k] for k, time in enumerate(row)] for machine, row in enumerate(idle_time)]

    # 计算总能耗
    return sum(sum(machine_energy)) + sum(sum(idle_energy)),

# 定义问题的类型（最小化两个目标）
creator.create("FitnessMulti", base.Fitness, weights=(-1.0, -1.0))
creator.create("Individual", list, fitness=creator.FitnessMulti)

# 定义遗传算法的参数
toolbox = base.Toolbox()
toolbox.register("attr_factory", random.randint, 0, f-1)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_factory, n)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("evaluate", lambda individual: (makespan(individual), energy(individual)))
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutUniformInt, low=0, up=f-1, indpb=0.1)
toolbox.register("select", tools.selNSGA2)

然后，使用NSGA-II算法求解该问题：

# 定义NSGA-II算法的参数
pop_size = 100  # 种群大小
ngen = 100  # 迭代次数
CXPB = 0.9  # 交叉概率
MUTPB = 0.1  # 变异概率

# 初始化种群
pop = toolbox.population(n=pop_size)

# 进化
for gen in range(ngen):
    # 生成子代
    offspring = tools.selTournamentDCD(pop, len(pop))
    offspring = [toolbox.clone(ind) for ind in offspring]
    for ind1, ind2 in zip(offspring[::2], offspring[1::2]):
        if random.random() < CXPB:
            toolbox.mate(ind1, ind2)
            del ind1.fitness.values
            del ind2.fitness.values
    for ind in offspring:
        if random.random() < MUTPB:
            toolbox.mutate(ind)
            del ind.fitness.values

    # 合并种群
    pop = tools.selBest(pop + offspring, pop_size)

# 输出最优解
best_ind = tools.selBest(pop, 1)[0]
print("Makespan = %.2f, Energy = %.2f" % (makespan(best_ind)[0], energy(best_ind)[0]))
print("Factory assignments:", best_ind)

该程序将输出求解结果，包括最小化完工时间和最小化总能耗的两个目标函数值，以及每个工件分配给哪个工厂。