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


  在橢圓x2+4y2=4 中有很多內接的矩形,這些矩形的邊平行於x軸和y軸,找出它們中面積最大的一個。

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

  矩形的兩條邊長分別是2x和2y,面積是A=4xy 。還知道xy都在橢圓上,因此問題可以轉換為約束條件下的極值:

數學方案

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

  由於拉格朗日乘子法無法確定最值的類型,所以還要計算函數邊界值。當P在橢圓上移動時,如果正好落在x軸上,則長方形退化成直線,此時面積是0;另一個邊界值是P落在y軸上,面積也是0。所以判定4是面積的最大值,P的坐標是(1.414, 0.707)

拋開數學的遺傳算法

  格朗日乘子法很簡單,但前提是需要了解偏導、梯度、拉格朗日乘子等一系列前置知識,而遺傳算法的優點是不需要復雜的數學知識,可以直接對問題編程求解,這對龐大的程序員群體來說無疑是個福音。

  可以通過橢圓得到xy的關系:

  矩形的面積可以寫成:

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

  由於xy都大於等於0,所以當4xy達到最大值時,(4xy)2也將達到最大。在求最值問題時,常系數起不到任何作用,因此求(4xy)2的最大值相當於求x2y2的最大值。結合xy的關系,原問題轉換為:

  現在看起來簡單多了。

  我們把問題精確到小數點后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位的”


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM