退而求其次(4)——橢圓中的最大矩形

　　在橢圓x²+4y²=4 中有不少內接的矩形，這些矩形的邊平行於x軸和y軸，找出它們中面積最大的一個。

　　先做圖，橢圓的中心在原點，其內接矩形的中心也在原點。設矩形的其中一點內接橢圓於P(x,y) ， P在第一象限：

　　矩形的兩條邊長分別是2x和2y，面積是A=4xy 。還知道x和y都在橢圓上，所以問題能夠轉換爲約束條件下的極值：

數學方案

　　最直接的方案是使用拉格朗日乘子法：

　　因爲拉格朗日乘子法沒法肯定最值的類型，因此還要計算函數邊界值。當P在橢圓上移動時，若是正好落在x軸上，則長方形退化成直線，此時面積是0；另外一個邊界值是P落在y軸上，面積也是0。因此斷定4是面積的最大值，P的座標是(1.414, 0.707)

拋開數學的遺傳算法

　　格朗日乘子法很簡單，但前提是須要了解偏導、梯度、拉格朗日乘子等一系列前置知識，而遺傳算法的優勢是不須要複雜的數學知識，能夠直接對問題編程求解，這對龐大的程序員羣體來講無疑是個福音。

　　能夠經過橢圓獲得x和y的關係：

　　矩形的面積能夠寫成：

　　只須要對x進行基因編碼就能夠利用遺傳算法求得最優解。在編寫代碼前，能夠考慮是否可以去掉影響算法運行速度和計算精度根號。

　　因爲x和y都大於等於0，因此當4xy達到最大值時，(4xy)²也將達到最大。在求最值問題時，常係數起不到任何做用，所以求(4xy)²的最大值至關於求x²y²的最大值。結合x和y的關係，原問題轉換爲：

　　如今看起來簡單多了。

　　咱們把問題精確到小數點後3位，從0.000到1.000之間有999個數，因爲999轉換成二進制是1111100111，所以可使用10位二進制基因編碼表示y。

　　代碼以下：

  1 import math
  2 import random
  3 import numpy as np
  4 import matplotlib.pyplot as plt
  5
  6 MAX, MIN = 0.999, 0 # y的最大值和最小值
  7 CODE_SIZE = 10 # 基因編碼長度
  8 POPULATION_SIZE = 20 # 種羣數量
  9
 10 def print_solution(code):
 11     '''
 12      打印解決方案
 13     :param code: 基因編碼
 14     '''
 15     x, y, A = decode(code)
 16     print('x = {0}, y = {1}, A = {2}'.format(x, y, A))
 17
 18 def decode(code:list):
 19     '''
 20     解碼
 21     :param code: 基因編碼
 22     :return: x,y,4xy （x,y都放大了1000倍）
 23     '''
 24     y = int(''.join(code), 2) * 0.001  # 將code轉換成十進制數
 25     if y > MAX:  # y超過了邊界
 26         return -1, -1, -1
 27     x = math.sqrt(4 - 4 * (y ** 2)) # x = sqrt(4 -4(y^2))
 28     # 使用退一法保留小數點後三位有效數字
 29     x, y = float('%.3f' % x), float('%.3f' % y)
 30     return x, y, 4 * x * y
 31
 32 def fitness_fun(code):
 33     ''' 適應度函數 '''
 34     return decode(code)[2]
 35
 36 def max_fitness(population):
 37     ''' 種羣中的最優適應個體 '''
 38     return max([fitness_fun(code) for code in population])
 39
 40 def init_population():
 41     ''' 構造初始種羣 '''
 42     population = []
 43     for i in range(POPULATION_SIZE):
 44         population.append(random.choices(['0', '1'], k=CODE_SIZE))
 45     return population
 46
 47 def selection(population):
 48     ''' 精英選擇策略'''
 49     # 按適應度從大到小排序
 50     pop_parents = sorted(population, key=lambda x: fitness_fun(x), reverse=True)
 51     # 選擇種羣中適應度最高的40%做爲精英
 52     return pop_parents[0: int(POPULATION_SIZE * 0.4)]
 53
 54 def crossover(population):
 55     ''' 單點交叉 '''
 56     pop_new = []  # 新種羣
 57     code_len = len(population[0])  # 基因編碼的長度
 58     for i in range(POPULATION_SIZE):
 59         p = random.randint(1, code_len - 1)  # 隨機選擇一個交叉點
 60         r_list = random.choices(population, k=2)  # 選擇兩個隨機的個體
 61         pop_new.append(r_list[0][0:p] + r_list[1][p:])
 62     return pop_new
 63
 64 def mutation(population):
 65     '''單點變異'''
 66     code_len = len(population[0])  # 基因編碼的長度
 67     mp = 0.2 # 變異率
 68     for i, r in enumerate(population):
 69         if random.random() < mp:
 70             p = random.randint(0, code_len - 1)  # 隨機變異點
 71             r[p] = '1' if r[p] == '0' else '1'
 72             population[i] = r
 73
 74 def ga():
 75     ''' 遺傳算法 '''
 76     population = init_population() # 構建初始化種羣
 77     max_fit = max_fitness(population) # 種羣最優個體的適應度
 78     max_fit_list = [max_fit]  # 每一代種羣的最優適應度
 79     i = 0
 80     while i < 5: # 若是連續5代沒有改進，結束算法
 81         pop_next = selection(population) # 選擇種羣
 82         pop_new = crossover(pop_next) # 交叉
 83         mutation(pop_new) # 變異
 84         max_fit_new = max_fitness(pop_new) # 新種羣中最優個體的適應度
 85         if max_fit < max_fit_new:
 86             max_fit = max_fit_new
 87             i = 0
 88         else:
 89             i += 1
 90         population = pop_new
 91         max_fit_list.append(max_fit_new)
 92     # 按適應度值從大到小排序
 93     population = sorted(population, key=lambda x: fitness_fun(x), reverse=True)
 94     # 返回最優的個體和每一代種羣中最優個體的適應度
 95     return population[0], max_fit_list
 96
 97 def pop_curve(max_fit_list):
 98     ''' 顯示種羣進化曲線 '''
 99     x = np.arange(1, len(max_fit_list) + 1, 1)
100     y = np.array(max_fit_list)
101     plot = plt.plot(x, y, '-')
102     plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
103     plt.title('種羣進化曲線')
104     plt.xlabel('種羣代數')
105     plt.ylabel('種羣總成本')
106     plt.show()
107
108 if __name__ == '__main__':
109     # code = ['1','0','1','1','0','0','0','0','1','1']  # 707的二進制基因編碼
110     # print_solution(code)
111     best, max_fit_list = ga()
112     print_solution(best)
113     pop_curve(max_fit_list)

　　decode方法先將y的基因編碼轉換成對應的十進制數，以後乘以0.001變成相應的小數。適應度函數的值是矩形的面積。遺傳算法的主體方法ga()除了獲得最優個體外，還額外返回了每一代種羣中最優個體的適應度，以便展現種羣進化曲線。在實際應用中，能夠經過觀察進化曲線來觀察算法在哪裏收斂，從而調整算法的終止條件。

　　種羣的進化曲線

　　能夠看到，算法收斂的至關快。一種可能的結果是：

　　

---

　　 做者：我是8位的

　　出處：http://www.cnblogs.com/bigmonkey

　　本文以學習、研究和分享爲主，如需轉載，請聯繫本人，標明做者和出處，非商業用途！ 

　　掃描二維碼關注公衆號「我是8位的」