遺傳算法010.解線性方程

2021-02-21 蝦神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<β上連續，則由上述通解公式可知
  • 人工智慧之遺傳算法(GA),搜索最優解的方法
    人工智慧之遺傳算法(GA),搜索最優解的方法 工程師8 發表於 2018-05-11 10:35:00 導讀：人們一提到遺傳算法（
  • 寫一個用迭代法解方程的Java程序
    歡迎點擊「算法與編程之美」關注我們！ 本文首發於微信公眾號："算法與編程之美"，歡迎關注，及時了解更多此系列文章。 2.解法介紹 牛頓迭代法是一種線性化方法，其基本思想是將非線性方程f(x)= 0逐步歸結-為某種線性方程來求解.設已知方程f(x)=0有近似根
  • 了解高階線性微分方程——初識二階線性微分方程
    小編照舊當大家都做了哦 ，現在微分方程篇已經算是複習了一半了，也不知道大家複習得怎麼樣，不過每天有進步終究還是好的。對於不想荒廢大學四年的同學，小編建議每天還是應該做一些學的課程裡面的題目，每天都有那份感覺在那裡，最終要用到的時候起碼不會生疏。下面小編開始對答案了。
  • 一階線性微分方程
    一階線性線性微分方程：形如的方程為一階線性微分方程。（1）當Q(x)≡0時，為一階線性齊次微分方程。（2）當Q(x)≠0時，則為一階線性非齊次微分方程。兩邊積分得 ln|y-3|=ln|x|+lnC即通解為y=Cx+3解 所給方程是一階線性齊次方程，且 P(x) = sin x，則
  • Python環境下的8種簡單線性回歸算法
    選自Medium作者：Tirthajyoti Sarkar機器之心編譯參與：晏奇、劉曉坤本文中，作者討論了 8 種在 Python 環境下進行簡單線性回歸計算的算法，不過沒有討論其性能的好壞，而是對比了其相對計算複雜度的度量。
  • 微分方程05 一階線性方程01
    叫做一階線性微分方程）方程。不同形式的微分方程解法不會像求導函數那樣具有有限固定的法則可依據，很多類微分方程問題都是一個相對獨立的孤島，歷經無數數學家的努力，很多微分方程或其中的某些特殊形式獲得了解析解， 而還有一些方程在現代計算機的幫助下獲得良好的數值解法。
  • 為什麼能用常數變易法求解一階非齊次線性微分方程
    記得在學高等數學微分方程這一章的時候，對於一階非齊次線性微分方程的計算，老師講過一種不用背公式的簡單算法——常數變易法。
  • 什麼是遺傳算法?怎樣繪製遺傳算法流程圖
    遺傳算法是一種通過模擬自然進化過程搜索最優解的操作方法，遺傳算法有三個基本算子選擇、交叉和變異。對於遺傳算法我們也可以使用流程圖對其整個過程進行總結歸納，那要怎樣繪製遺傳算法流程圖呢？下面是分享的簡單操作方法，希望可以幫助大家。
  • 玩轉線性代數(21)線性方程組解的判斷與求法
    討論內容：線性方程組的解的判斷與矩陣的秩，矩陣方程的解法記錄：7.1線性方程組解的討論我：「上節課真夠費腦子的，今天好多了，引入了矩陣的秩後，下面來討論一下矩陣的秩與線性方程組的解的判斷．」小明：「如何判斷矩陣方程Ax=B是否有解呢？」我：「這個問題好，我來分析一下．」問題：AX=B的解的存在性．
  • 2016考研數學之解線性方程組的方法
    線性方程組的三種形式包括原始形式、矩陣形式、向量形式，高斯消元法是最基礎和最直接的求解線性方程組的方法，其中涉及到三種對方程的同解變換：（1）把某個方程的k倍加到另外一個方程上去；（2）交換某兩個方程的位置；（3）用某個常數k乘以某個方程。我們把這三種變換統稱為線性方程組的初等變換。
  • 詳解遺傳算法
    遺傳算法是受達爾文的進化論的啟發，借鑑生物進化過程而提出的一種啟發式搜索算法。因此在介紹遺傳算法前有必要簡單的介紹生物進化知識。 遺傳算法思想 　　借鑑生物進化論，遺傳算法將要解決的問題模擬成一個生物進化的過程，通過複製、交叉、突變等操作產生下一代的解，並逐步淘汰掉適應度函數值低的解，增加適應度函數值高的解。這樣進化N代後就很有可能會進化出適應度函數值很高的個體。
  • 系統微分方程的解―系統的全響應
    一、 線性系統微分方程線性的證明本文引用地址：http://www.eepw.com.cn/article/155509.htm線性系統必須同時滿足齊次性與疊加性所以，要證明線性系統的微分方程是否是線性的，就必須證明它是否同時滿足齊次性與疊加性。
  • MATLAB入門之MATLAB求線性方程組的解
    求解多元一次線性方程組的解，可以轉化為求解矩陣解的問題。例如：求下列線性方程組的解解：此方程可列成兩組不同的矩陣方程形式。一是，設X=[x1;x2;x3;x4]為列向量，矩陣A= [1 4 –7 6;0 2 1 1;0 1 1 3;1 0 1 –1]，B=[0;-8;-2;1]為列向量，則方程形式為AX=B，其求解過程用左除：由此可見，A\B 的確與inv(A)*B 相等。
  • 【物理數學】解微分方程為什麼會出現個e?
    Q：解微分方程為什麼會出現個e？
  • 多元方程的解如何表示?
    說到方程，最簡單的莫過於一元一次方程：只有一個等式，又沒有複雜運算，只靠基本線性運算就可以算出答案。一元一次方程的名字中包含元和次，元是指未知數的個數，而次是指未知數最高有幾次冪。既然有一元一次方程，就當然有多元一次方程，一元多次方程和多元多次方程。多元多次方程對應的函數其實是冪函數，並未涉及到其他複雜運算，更普遍的形式是多元方程，方程裡面包含各種運算（指數運算，對數運算，三角函數、反三角函數，甚至微分和積分運算等等）多元方程的解很可能有無數個解，其實多元方程的解只是對應多元函數的自變量值。
  • 2017考研數學之微分方程解的結構
    微分方程是高等數學的重要組成部門，無論對於數一數二，還是數三，每一年都會考察，並且題型不定，會出現選擇題，填空題，也會出現解答題。對於微分方程這一部分，總的解題思路是首先判斷微分方程的形式和結構，如是否為線性方程，齊次方程等等，是一階、二階還是更高階的方程；然後根據方程的形式和結構來求解方程。
  • 常微分方程的級數解
    在科學研究中，經常需要解含有未知函數的導數的方程，這就是微分方程。