遺傳算法010.解線性方程

2021-02-19 蝦神daxialu
寫在前面的話

解方程在實際應用中是非常常見的問題——相信每個同學在讀書的時候都給數學方程給折磨過。與數學領域的解方程不同，數學領域可以直接用分數或者無理數來表達，比如：

但是在計算機領域，就必須表達為：

import math
a = math.sqrt(2)
b = 2/3
print("a = {}  b = {} a*b = {}".format(a,b,a*b))
# a = 1.4142135623730951  b = 0.6666666666666666 a*b = 0.9428090415820634
從上面可以看見，計算機對於循環小數或者無理數，採取的是精度截取的方式來實現，所以計算機裡面得到的永遠是近似值。
回到解方程裡面，與上面的計算一樣，數學裡面的解方程可以得到表達式組合，而計算機裡面最終會給出近似值來，所以計算機解方程與（中學）數學原理是不太一樣的。
對於解方程的方法，一般有兩種：
第一種稱之為數學方法，其中包括了：
LU分解(LU Factorization)
高斯消元法
QR（正交三角）分解法
第二種就是所謂迭代法，也就是計算機常用的方法，而我們要介紹的遺傳算法解方程，也屬於迭代法中的一種。
在本節中，我們將看到，遺傳算法如何來解線性方程組。
比如我們現在有這樣一組方程：
x + 2y = 4
4x + 4y = 12
用高斯消元法，很容易就能算出來：
方程2 - （方程1*2）得：2x = 4 解得x = 2 帶入得 y = 1
同樣的Python裡面也有包括numpy和scipy這樣的專業計算工具包，可以來用解方程：
import numpy
A = numpy.array([[1,2],[4,4]])
B = numpy.array([4, 12])
inv = numpy.linalg.inv(A).dot(B)
print("x1 = {}  x2 = {}".format(inv[0],inv[1]))
# x1 = 2.0  x2 = 1.0
算法實現
前面是通用方法，不解釋了，大家看以前的文章：本節用了一個Python的內部包fration（分數包），這個包可以實現直接保留為分數，而不會計算成浮點數，比如：
# 三分之一乘以10
# 傳統情況下，會直接轉化為浮點數：
print(1/3 * 10.0)
# 3.333333333333333

#用fration包，可以保留分數計算
print(Fraction(1,3) * 10)
# 10/3
import random,datetime
from bisect import bisect_left
from math import exp
from fractions import Fraction

# 初始化初代種群的方法
def _generate_parent(length, geneSet, get_fitness):
    genes = []
    while len(genes) < length:
        sampleSize = min(length - len(genes), len(geneSet))
        genes.extend(random.sample(geneSet, sampleSize))
    fitness = get_fitness(genes)
    return Chromosome(genes, fitness)

# 封裝的進化方法
def _mutate_custom(parent, custom_mutate, get_fitness):
    childGenes = parent.Genes[:]
    custom_mutate(childGenes)
    fitness = get_fitness(childGenes)
    return Chromosome(childGenes, fitness)

#進化執行
def get_best(get_fitness, targetLen, optimalFitness, geneSet, display,
             custom_mutate=None, custom_create=None, maxAge=None):
    if custom_mutate is None:
        def fnMutate(parent):
            return _mutate(parent, geneSet, get_fitness)
    else:
        def fnMutate(parent):
            return _mutate_custom(parent, custom_mutate, get_fitness)

    if custom_create is None:
        def fnGenerateParent():
            return _generate_parent(targetLen, geneSet, get_fitness)
    else:
        def fnGenerateParent():
            genes = custom_create()
            return Chromosome(genes, get_fitness(genes))

    for improvement in _get_improvement(fnMutate, fnGenerateParent, maxAge):
        display(improvement)
        if not optimalFitness > improvement.Fitness:
            return improvement


def _get_improvement(new_child, generate_parent, maxAge):
    parent = bestParent = generate_parent()
    yield bestParent
    historicalFitnesses = [bestParent.Fitness]
    while True:
        child = new_child(parent)
        if parent.Fitness > child.Fitness:
            if maxAge is None:
                continue
            parent.Age += 1
            if maxAge > parent.Age:
                continue
            index = bisect_left(historicalFitnesses, child.Fitness, 0,
                                len(historicalFitnesses))
            proportionSimilar = index / len(historicalFitnesses)
            if random.random() < exp(-proportionSimilar):
                parent = child
                continue
            bestParent.Age = 0
            parent = bestParent
            continue
        if not child.Fitness > parent.Fitness:
            # same fitness
            child.Age = parent.Age + 1
            parent = child
            continue
        child.Age = 0
        parent = child
        if child.Fitness > bestParent.Fitness:
            bestParent = child
            yield bestParent
            historicalFitnesses.append(bestParent.Fitness)

#種群類
class Chromosome:
    def __init__(self, genes, fitness):
        self.Genes = genes
        self.Fitness = fitness
        self.Age = 0
下面是健壯性驗證的方法：
把每個迭代出來的未知數的值，帶入到方程裡面，如果等於0，這表示通過：
把x + 2y = 4 寫成 x + 2y - 4 = 0這樣的模式，如果帶入未知數，最後等於0則表示通過。
def get_fitness(genes, equations):
    fitness = Fitness(sum(abs(e(genes)) for e in equations))
    return fitness
顯示方法，不解釋：
def display(candidate, startTime, fnGenesToInputs):
    timeDiff = datetime.datetime.now() - startTime
    symbols = "xyza"
    result = ', '.join("{} = {}".format(s, v)
                       for s, v in
                       zip(symbols, fnGenesToInputs(candidate.Genes)))
    print("{}\t{}\t{}".format(
        result,
        candidate.Fitness,
        timeDiff))
進化方法：首先選擇一個基因進行變化，然後通過改變窗口的大小的方法，不斷嘗試進行基因替換。每次通過循環，從選擇的索引中提取下一個基因，並使用窗口大小計算出需要搜索範圍。這樣，就可以限制可以挑選的基因數量，如果到達數組的一端，就停止。最後，用從新的邊界範圍內隨機抽取的一個基因取代當前的基因。
def mutate(genes, sortedGeneset, window, geneIndexes):
    indexes = random.sample(geneIndexes, random.randint(1, len(genes))) \
        if random.randint(0, 10) == 0 else [random.choice(geneIndexes)]
    window.slide()
    while len(indexes) > 0:
        index = indexes.pop()
        genesetIndex = sortedGeneset.index(genes[index])
        start = max(0, genesetIndex - window.Size)
        stop = min(len(sortedGeneset) - 1, genesetIndex + window.Size)
        genesetIndex = random.randint(start, stop)
        genes[index] = sortedGeneset[genesetIndex]
使用窗口，主要是為了讓我們在尋找數據的時候，儘量接近目前要替換的基因，比如要替換的是個正數，那麼替換它的也儘可能是一個正數，反之亦然。
這樣我們就可以慢慢改變當前最佳基因周圍的搜索範圍。這使得算法能夠比較容易的找到下一個最優的改進可能。
兩個常規類：健壯性類和窗口類：
class Fitness:
    def __init__(self, totalDifference):
        self.TotalDifference = totalDifference

    def __gt__(self, other):
        return self.TotalDifference < other.TotalDifference

    def __str__(self):
        return "diff: {:0.2f}".format(float(self.TotalDifference))


class Window:
    def __init__(self, minimum, maximum, size):
        self.Min = minimum
        self.Max = maximum
        self.Size = size

    def slide(self):
        self.Size = self.Size - 1 if self.Size > self.Min else self.Max
下面就是解題的過程了，一個封裝了一堆內部方法的方法：
def solve_unknowns(numUnknowns, geneset, equations,fnGenesToInputs):
        startTime = datetime.datetime.now()
        maxAge = 50
        window = Window(max(1, int(len(geneset) / (2 * maxAge))),
                        max(1, int(len(geneset) / 3)),
                        int(len(geneset) / 2))
        geneIndexes = [i for i in range(numUnknowns)]
        sortedGeneset = sorted(geneset)

        def fnDisplay(candidate):
            display(candidate, startTime, fnGenesToInputs)

        def fnGetFitness(genes):
            return get_fitness(genes, equations)

        def fnMutate(genes):
            mutate(genes, sortedGeneset, window, geneIndexes)

        optimalFitness = Fitness(0)
        best = get_best(fnGetFitness, numUnknowns, optimalFitness,
                                geneset, fnDisplay, fnMutate, maxAge=maxAge)
然後是接二元一次方程：
用遺傳算法解方程，最大的難題，就是在於如何構建基因庫。比如我們這裡做二元一次方程，設計的基因就是從-5到5這10個整數。如果你的解不在基因庫裡面，那是肯定解不出來的。
def test_2_unknowns():
    geneset = [i for i in range(-5, 5) if i != 0]

    def fnGenesToInputs(genes):
        return genes[0], genes[1]

    def e1(genes):
        x, y = fnGenesToInputs(genes)
        return x + 2 * y - 4

    def e2(genes):
        x, y = fnGenesToInputs(genes)
        return 4 * x + 4 * y - 12
    equations = [e1, e2]
    solve_unknowns(2, geneset, equations, fnGenesToInputs)
運行結果如下：
x = 1, y = -5	diff: 41.00	0:00:00
x = 3, y = -3	diff: 19.00	0:00:00
x = 4, y = -3	diff: 14.00	0:00:00
x = 4, y = -1	diff: 2.00	0:00:00.000999
x = 1, y = 2	diff: 1.00	0:00:00.002992
x = 2, y = 1	diff: 0.00	0:00:00.018948
同樣，我們可以簡單的改一改，比如把方程改成：
30x + 15y = 600
15x + 20y = 200
那麼初始化解題的方法就要改成：
首先我們不確定基因庫範圍的時候，給一個比較高的（這裡還是已經知道接肯定是整數，否則也不行） 之後改兩個方程的構建。
def test_2_unknowns2():
    geneset = [i for i in range(-100,100) if i != 0]

    def fnGenesToInputs(genes):
        return genes[0], genes[1]

    def e1(genes):
        x, y = fnGenesToInputs(genes)
        return 30 *x + 15 * y - 600

    def e2(genes):
        x, y = fnGenesToInputs(genes)
        return 15 * x + 20 * y - 200
    equations = [e1, e2]
    solve_unknowns(2, geneset, equations, fnGenesToInputs)
test_2_unknowns2()
運行結果：
x = 22, y = -26	diff: 720.00	0:00:00
x = 27, y = -26	diff: 495.00	0:00:00.001004
x = 45, y = -26	diff: 405.00	0:00:00.001004
x = 42, y = -26	diff: 360.00	0:00:00.001004
x = 41, y = -26	diff: 345.00	0:00:00.001004
x = 32, y = -26	diff: 270.00	0:00:00.001004
x = 14, y = 5	diff: 215.00	0:00:00.001004
x = 14, y = -1	diff: 205.00	0:00:00.001993
x = 20, y = -1	diff: 95.00	0:00:00.001993
x = 20, y = -2	diff: 90.00	0:00:00.001993
x = 20, y = -3	diff: 85.00	0:00:00.002998
x = 21, y = -3	diff: 70.00	0:00:00.002998
x = 21, y = -5	diff: 60.00	0:00:00.002998
x = 22, y = -5	diff: 45.00	0:00:00.002998
x = 23, y = -6	diff: 25.00	0:00:00.009972
x = 23, y = -7	diff: 20.00	0:00:00.018961
x = 24, y = -8	diff: 0.00	0:00:00.135640
下面這個帶有分數三元一次方程：
6x - 2y + 8z = 20
y + 8x * z = -1
2z * 6/x + 3y/2 = 6
構造初始化方法如下：
構建的解是從-5/4 到 5/4
不要問我如果真實未知，不可能知道範圍怎麼辦……這裡的作用是給大家說明算法，不是真實使用的方法，千萬別再真實工作中去使用這種方法來解方程，沒意義
def test_3_unknowns():
        geneRange = [i for i in range(-5, 5) if i != 0]
        geneset = [i for i in set(
            Fraction(d, e)
            for d in geneRange
            for e in geneRange if e != 0)]

        def fnGenesToInputs(genes):
            return genes

        def e1(genes):
            x, y, z = genes
            return 6 * x - 2 * y + 8 * z - 20

        def e2(genes):
            x, y, z = genes
            return y + 8 * x * z + 1

        def e3(genes):
            x, y, z = genes
            return 2 * z * Fraction(6, x) \
                   + 3 * Fraction(y, 2) - 6

        equations = [e1, e2, e3]
        solve_unknowns(3, geneset, equations, fnGenesToInputs)
執行結果如下：
x = 2/5, y = -4, z = -3/4	diff: 55.50	0:00:00
x = 2/5, y = -3, z = -3/4	diff: 55.00	0:00:00.000995
x = -3/5, y = -3, z = -3/4	diff: 29.70	0:00:00.000995
x = -3/5, y = -5, z = -3/4	diff: 21.50	0:00:00.001963
x = -2/3, y = -5, z = -3/4	diff: 20.00	0:00:00.005951
x = -3/5, y = -5, z = -2/3	diff: 19.90	0:00:00.019914
x = -1/4, y = -4, z = -1/4	diff: 18.00	0:00:00.036870
x = -1/5, y = -4, z = -1/5	diff: 17.48	0:00:00.045844
x = 1/2, y = -4/5, z = 1/5	diff: 17.20	0:00:00.059808
x = 1/3, y = -4/5, z = 1/5	diff: 15.53	0:00:00.059808
x = 1/3, y = -5/4, z = 1/5	diff: 14.86	0:00:00.061280
x = 1/3, y = -3/2, z = 1/5	diff: 14.48	0:00:00.063252
x = 5/3, y = -3/2, z = 1/5	diff: 14.38	0:00:00.067241
x = 3/4, y = -5, z = 3/4	diff: 2.50	0:00:00.068236
x = 3/5, y = -5, z = 3/4	diff: 2.30	0:00:00.071230
x = 3/5, y = -5, z = 2/3	diff: 2.03	0:00:00.076215
x = 2/3, y = -5, z = 3/4	diff: 0.00	0:00:00.092175
當然， 你還可以改成四元甚至五元方程，關鍵點在於你要能夠構建合理的基因庫，否則遺傳算法就沒有辦法使用了。
總結：
遺傳算法需要構建良好的基因庫，也需要有很好的終止條件，如果這兩個前提做不到的話，遺傳算法就沒辦法使用。

相關焦點

  • 考研數學:線性微分方程解的線性組合分析
    微分方程是考研數學的一個必考點，並且還常常出兩道題，佔十幾分，因此大家一定重視這一章，掌握好其基本知識和解題方法。微分方程中最重要的是線性微分方程，主要包括一階和二階線性微分方程，而線性微分方程又可分為線性齊次和非齊次線性微分方程，它們的解的線性組合是否仍為其方程的解?
  • 拆分——線性微分方程的解的結構
    話不多說，這篇文章算是微分方程這一章難點的開頭了。我們現在來複習線性微分方程的解的結構。這裡主要討論二階線性方程，並且考試中也只會出現二階，不會考到三階及其以上線性方程的。接下來先講齊次方程的解的結構：y''+P(x)y'+Q(x)y=0補充：線性相關及無關就是說如果所有的k都為0，那麼就是線性無關，否則就是線性相關。
  • 誤差的分析與減少及Matlab解線性方程的四種方法
    1、誤差的來源     模型誤差：數學模型與實際問題之間的誤差     觀測誤差：測量數據與實際數據的誤差     方法誤差：數學模型的精確解與數值方法得到的數值解之間的誤差
  • 常微分方程:線性微分方程解的三個重要特徵
    前一篇《帶你走進微積分的堂學習：一階線性微分方程式的基礎原理》詳細討論了線性微分方程的結構以及通解特性，本篇我們藉此機會指出一階線性微分方程解的三個重要特徵1）有一階線性微分方程>的通解是可以看出，它等於（1）的一個特解（對應於上式的C=0）再加相應的齊次線性（2）的通解，因此如果求得非齊次線性微分方程（1）的一個特解為y=φ1(x)和相應的齊次線性方程（2）的通解，則(1)的通解為2）設a(x)和b(x)在區間α<x<β上連續，則由上述通解公式可知
  • 線性代數拾遺(一):線性方程組、向量方程和矩陣方程
    線性方程組、向量方程和矩陣方程一、線性方程組線性代數，最基本的問題，就是解線性方程組了。線性方程組就是一組形如 a1x1+a2x2+⋯+anxn=ba1x1+a2x2+⋯+anxn=b的方程。比如，將增廣矩陣的第一、二行對換，那麼同時，它所代表的線性方程組中，第一、二個方程也進行了對調。1.2 線性方程組的解解一個線性方程組，就是通過對其矩陣形式行變換（三種方式：交換方程的先後順序，一個方程左右同乘以某數，和兩個方程相加） 轉換為行階梯形式。
  • 2017考研數學:n階線性微分方程的通解公式分析
    關於線性微分方程的通解公式，在一般高等數學教材中只是簡單地做了些介紹，並沒有進行詳細的分析證明，因此有很多同學對其感到有些困惑，對其含義和作用也不能很好理解，為了幫助2017考研學子消除這些困惑，本文對n階線性微分方程的通解公式做些分析和證明，供同學參考。
  • 2017考研數學:討論關於線性方程組的解的問題
    線性方程組是線性代數的考研數學線性代數重難點，也是考察重點，尤其是求解的問題，下面新東方網考研頻道和大家一起來討論下。 　　線性方程組的特點：方程是未知數的一次齊次式，方程組的數目s和未知數的個數n可以相同，也可以不同。
  • 今日PNAS:利用新型存儲器陣列一步解線性方程組和特徵向量
    線性代數中最核心的一個操作是解不同的矩陣方程，包括解線性方程組和特徵向量等。在傳統計算機上，解矩陣方程需要一些精心設計的算法，如高斯消元法、LU分解法，然後在多項式時間內（比如O(N3)，N是矩陣行/列數）獲得方程解。
  • 了解高階線性微分方程——初識二階線性微分方程
    小編照舊當大家都做了哦 ，現在微分方程篇已經算是複習了一半了，也不知道大家複習得怎麼樣，不過每天有進步終究還是好的。對於不想荒廢大學四年的同學，小編建議每天還是應該做一些學的課程裡面的題目，每天都有那份感覺在那裡，最終要用到的時候起碼不會生疏。下面小編開始對答案了。
  • 微分方程06 一階線性方程02
    一階線性微分方程有標準形式：
  • 一階線性微分方程
    的方程為一階線性微分方程。（1）當Q(x)≡0時，為一階線性齊次微分方程。（2）當Q(x)≠0時，則為一階線性非齊次微分方程。求解微分方程的通解解：齊次方程，分離變量得兩邊積分得 ln|y-3|=ln|x|+lnC即通解為y=Cx+3解 所給方程是一階線性齊次方程，且 P(x) = sin x，則
  • 2017考研線代重點:高斯消元法解線性方程組
    線性方程組的三種形式包括原始形式、矩陣形式、向量形式，高斯消元法是最基礎和最直接的求解線性方程組的方法，其中涉及到三種對方程的同解變換： 　　(1)把某個方程的k倍加到另外一個方程上去; 　　(2)交換某兩個方程的位置; 　　(3)用某個常數k乘以某個方程。我們把這三種變換統稱為線性方程組的初等變換。
  • 微分方程05 一階線性方程01
    叫做一階線性微分方程）方程。不同形式的微分方程解法不會像求導函數那樣具有有限固定的法則可依據，很多類微分方程問題都是一個相對獨立的孤島，歷經無數數學家的努力，很多微分方程或其中的某些特殊形式獲得了解析解， 而還有一些方程在現代計算機的幫助下獲得良好的數值解法。
  • 遺傳算法簡介、基本原理及算法實現
    遺傳算法算是一個比較複雜的算法，本文主要分為三個部分，第一部分遺傳算法的發展、功能、主要思想以及簡單的代碼命令；第二部分詳細介紹遺傳算法的機理，第三部分舉例並給出一個簡單地使用方法。雖然遺傳算法存在諸如陷入局部最優解，收斂速度緩慢等問題，人們也進行了很多的修改，但是鑑於遺傳算法本身就具有很高的性能，而且各種修正方案都存在一定的複雜性和非普適性，所以經典遺傳算法在應用領域還是幾乎絕對的主流。二、遺傳算法的原理遺傳的算法的背景是達爾文的進化論，沒錯！
  • Python環境下的8種簡單線性回歸算法
    選自Medium作者：Tirthajyoti Sarkar機器之心編譯參與：晏奇、劉曉坤本文中，作者討論了 8 種在 Python 環境下進行簡單線性回歸計算的算法有關各種線性回歸方法的代碼可以參閱筆者的 GitHub。其中大部分都基於 SciPy 包SciPy 基於 Numpy 建立，集合了數學算法與方便易用的函數。通過為用戶提供高級命令，以及用於操作和可視化數據的類，SciPy 顯著增強了 Python 的交互式會話。
  • 考研數學:線性方程組與線性微分方程的通解對比
    線性方程組是線性代數中的一個重要知識點，而線性微分方程是高等數學中微分方程部分的一個重要知識點，二者雖然分別屬於不同的數學課程內容，但其通解形式卻有著驚人的相似之處，有些同學在學習中感覺到了二者有相似之處，但並不十分清楚其相似在何處和怎麼相似，以及線性微分方程的通解是否包含其全部解，對此的蔡老師就這些問題做些歸納總結，供複習2018考研數學的同學和學習高等數學及線性代數的同學參考
  • 2017考研高數六大基本題型:微分方程解常微分方程
    新東方網>大學教育>考研>複習指導>數學>正文2017考研高數六大基本題型：微分方程解常微分方程 2016-10-26 16:48 來源：新東方網整理
  • 2017考研數學:二階常係數線性齊次差分方程的通解分析
    新東方網>大學教育>考研>複習指導>數學>正文2017考研數學：二階常係數線性齊次差分方程的通解分析 2016-03-30 09:37 來源：文都考研
  • (數一)二階常係數線性微分方程
    正文：什麼是二階常係數線性微分方程？二階常係數線性微分方程是形如y''+py'+qy=f(x)的微分方程，其中p，q是實常數。自由項f(x)為定義在區間I上的連續函數，即y''+py'+qy=0時，稱為二階常係數齊次線性微分方程。
  • 什麼是遺傳算法?怎樣繪製遺傳算法流程圖
    遺傳算法是一種通過模擬自然進化過程搜索最優解的操作方法，遺傳算法有三個基本算子選擇、交叉和變異。對於遺傳算法我們也可以使用流程圖對其整個過程進行總結歸納，那要怎樣繪製遺傳算法流程圖呢？下面是分享的簡單操作方法，希望可以幫助大家。