# !/usr/bin/env python
# -*- coding: utf-8 -*-

#  Reference: https://zhuanlan.zhihu.com/p/127527177
#  Usage: Example for solving multi-objective problem with geatpy.

import numpy as np
import geatpy as ea

class MyProblem(ea.Problem):

       # --------- 定义参数 --------- #
       def __init__(self):
              name ='BNH'
              M = 2 # 定义目标维数
              maxormins = [1] * M
              Dim = 2 # 定义决策变量维数
              varTypes = [0] * Dim # 定义决策变量的类型，其中0表示实数，1表示整数
              lb = [0] * Dim # 定义决策变量下界
              ub = [5, 3] # 定义决策变量上界
              lbin = [1] * Dim # 定义决策变量下边界
              ubin = [1] * Dim #决策变量上边界
              ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)

       # --------- 定义目标函数和约束 --------- #
       def aimFunc(self, pop):
              Vars = pop.Phen
              x1 = Vars[:, [0]]  # 定义决策变量矩阵
              x2 = Vars[:, [1]]
              f1 = 4 * x1 ** 2 + 4 * x2 ** 2
              f2 = 4 * (x1 - 5) ** 2 + 4 * (x2 - 5) ** 2
              pop.CV = np.hstack(
              [(x1 - 5)**2 + x2**2 - 25,
              -(x1 - 8)**2 - (x2 - 3)**2 + 7.7])  # 采用可行性法则处理约束
              pop.ObjV = np.hstack([f1, f2]) # 把求得的目标函数值赋值给种群pop的ObjV

       # --------- 计算全局最优解 --------- #
       def calReferObjV(self):
              N = 10000
              x1 = np.linspace(0, 5, N)
              x2 = x1.copy()
              x2[x1 >= 3] = 3
              return np.vstack((4 * x1 ** 2 + 4 * x2 ** 2, 4 * (x1 - 5) ** 2 + 4 * (x2 - 5) **2)).T

# --------- 实例化问题对象 --------- #
problem = MyProblem()

# --------- 种群设置 --------- #
Encoding ='RI' # 定义编码方式
NIND = 100 # 定义种群规模
Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges, problem.borders) # 创建区域描述器
population = ea.Population(Encoding, Field, NIND) # 实例化种群对象

# --------- 参数设置 --------- #
myAlgorithm = ea.moea_NSGA2_templet(problem, population) # NSGA2是一个多目标优化算法
myAlgorithm.MAXGEN = 200 # 定义最大遗传代数

# --------- 模型优化 --------- #
NDSet = myAlgorithm.run()
print('用时：%f秒'%(myAlgorithm.passTime))
#计算指标
PF = problem.getReferObjV()
if PF is not None and len(NDSet) != 0:
       IGD = ea.indicator.IGD(NDSet[1].ObjV, PF) # 计算IGD(Inverted Generational Distance,反转世代距离)用于评价算法收敛性和多样性
       print('IGD:%f'%IGD)