线性与非线性规划：复合形方法

S_Radar

已于 2022-09-20 11:23:09 修改

阅读量1.5k

点赞数 1

分类专栏：最优化线性与非线性规划文章标签： python

于 2022-08-12 20:29:38 首次发布

本文链接：https://blog.csdn.net/qq_45397810/article/details/126309240

版权

最优化同时被 2 个专栏收录

10 篇文章 1 订阅

订阅专栏

线性与非线性规划

10 篇文章 1 订阅

订阅专栏

文章目录

前言
一、复合形方法
二、算法实现
总结

前言

本系列涉及线性与非线性规划中的几种规划算法

1.本节介绍复合形方法（单纯形搜索方法）
2.参考：陈宝林-最优化理论与算法
3.采用python编程实现，已测试，代码可行

一、复合形方法

这里的复合形方法，我觉得和参考书中的单纯形搜索方法是一个意思。基本思想是先生成k个可行点，得到一个初始的复合形，然后改变复合形的形状（反射、扩张、收缩、压缩）。
这考虑反射这一种手段的复合形方法的具体细节参见参考书和下述。

二、算法实现

1.算法流程

只是反射方式下的复合形方法的流程如下图：
算法步骤
算法流程框图

2.算法细节

初始复合形

反射机理
反射点位置
扩张机理
收缩机理
压缩机理

3.代码

#复合型方法，随机初始区间X0_R=[[a,b],[],[]...]未给出，Xi自变量，X0给定一个初始可行点初值
#定点数n+1<=k<=2n,收敛精度delta,a0每次缩减的倍数，lamda,a0的最小限制a0_limit=e-10,反射a0=1.3，b0，扩张a1=1,b1,缩小a2=0.7,b2,压缩a3=0.5
import math
from sympy import *
import random
import numpy as np
class Complex_method:
    def __init__(self, F, G, Xi, X0_R, k, X0=None):
        self.f = F
        self.G = G
        self.X0_R = X0_R
        self.i = len(X0_R)# 自变量个数
        self.Xi = Xi  # 定义自变量
        self.k = k  # 产生的顶点数
        if not X0:
            self.Xk = np.array(self.random_X(self.k,[]))
            #self.Xk = np.array([[7,8],[6,8],[8,7]],dtype=float)
        else:
            Xk = []
            Xk.append(X0)
            self.Xk = np.array(self.random_X(k-1,Xk))
    # 生成自变量字典，如何实现字典的numpy向量计算
    def x_dict(self,X):
        x = {}
        for index,item in enumerate(self.Xi):
            x[item] = X[index]
        self.X_dict = x
    #随机顶点Xk是存放顶点的集合
    def random_X(self,k,Xk):
        n = 0
        while n < k:
            X = []
            for index,item in enumerate(self.Xi):
                X.append(random.uniform(self.X0_R[index][0],self.X0_R[index][1]))
            if self.condition_G(X):
                Xk.append(X)
                n += 1
        return Xk
    #约束条件
    def condition_G(self,X):
        self.x_dict(X)
        n = True
        for i in self.G:
            if i.subs(self.X_dict) > 0:
                n = false
        return n
    #计算目标函数值
    def f_value(self,X):
        self.x_dict(X)
        self.fX = self.f.subs(self.X_dict)
        return self.fX
    #求各顶点的函数值并排序
    def fXk_value(self):
        fXk = []
        for i in self.Xk:
            fXk.append(self.f_value(i))
        self.fXk = np.array(fXk)#顶点的值
        self.fmeans = self.fXk.mean()#平均值
        max_fxk = max(fXk)
        min_fxk = min(fXk)
        acc = math.sqrt(sum((np.array(fXk)-min_fxk)**2)/(self.k-1))#
        self.fXl = min_fxk
        self.fXH = max_fxk
        max_index = fXk.index(max_fxk)
        min_index = fXk.index(min_fxk)
        self.XL = self.Xk[min_index]#最小(好)
        self.XL_index = min_index
        self.XH = self.Xk[max_index]#最大(坏)
        self.XH_index = max_index
        del fXk[max_index]
        max_1 = max(fXk)
        self.fXG = max_1
        max1_index = self.fXk.tolist().index(max_1)
        self.XG = self.Xk[max1_index]#次大(坏)
        self.XG_index = max1_index
        return acc
    #计算中心点和
    def center(self):
        Xcenter = (np.sum(self.Xk,axis=0) - self.XH) /(self.k-1)
        return Xcenter
    #计算反射点
    def reflect(self,XC,a0):
        XR = XC + a0 * (XC-self.XH)
        #self.Xk[self.XH_index] = XR
        return XR
    #扩张a1
    def again_reflect(self,XR,XC,a1):
        n = false
        XE = XR + a1 * (XR-XC)
        if self.condition_G(XE):  # 扩张后的点可行
            if self.f_value(XE) < self.f_value(XR):
                self.Xk[self.XH_index] = XE
                n = True
        return  n
    #缩小
    def reduct(self,XC,a2):
        n = false
        XK = self.XH + a2 * (XC - self.XH)
        if self.condition_G(XK):  # 扩张后的点可行
            if self.f_value(XK) < self.f_value(self.XH):
                self.Xk[self.XH_index] = XK
                n = True
        return n
    #压缩
    def compress(self,XL,Xk,a3):
        Xk = XL - a3 * (XL - Xk)
        Xk[self.XL_index] = XL
        self.Xk = Xk
        return Xk
    #开始计算
    #收敛精度delta, 每次缩减的倍数lamda,反射a0 = 1.3, a0的最小限制a0_limit = e - 10，b0，扩张a1 = 1, b1, 缩小a2 = 0.7, b2, 压缩a3 = 0.5
    def find_(self,delta,lamda,a0,a0_limit):
        a1 = a0
        aa = 0
        bb = 0
        ci = 0
        while True:
            acc = self.fXk_value()#是否为最优解
            print(f'顶点为Xk={self.Xk}')
            print(f'XL={self.XL},fXL={self.fXl},XG={self.XG},fXG={self.fXG},XH={self.XH},fXH={self.fXH}')
            print(f'精度acc={acc}')
            if  acc<= delta:
                print(f'最优解为x*={self.XL},极小值为f*={self.fXl}')
                break
            else:
                Xcenter = self.center()  # 中心可行否？
                a0 = a1#每次重新反射都要还原a0
                while True:
                    if aa == 1:#为了
                        aa = 0  # 退出到acc
                        break#即用来控制elif self.condition_G(Xcenter)的退出
                    elif self.condition_G(Xcenter) :
                        print(f'可行中心点={Xcenter}')
                        while True:
                            XR = self.reflect(Xcenter,a0)#反射点行否？
                            if self.condition_G(XR) :
                                while True:
                                    bb = 0
                                    if self.f_value(XR) < self.f_value(self.XH):#是否下降
                                        ci += 1
                                        self.Xk[self.XH_index] = XR
                                        print(f'反射点XR={XR}')
                                        aa = 1
                                        break#控制从if self.f_value(XR) < self.f_value(self.XH):的退出

                                    else:#不下降
                                        if a0 <= a0_limit:
                                            self.XH = self.XG###关于XG是否也要换为XH？？交换的效果更好
                                            print(f'a0达到极限')
                                            break#控制从if a0 <= a0_limit:退出
                                        else:
                                            a0 *= lamda#
                                            print(f'不下降的反射点XR={XR},a0减小，a0={a0}')
                                            bb = 1
                                            break
                                if bb == 1:
                                    continue#控制从a0 *= lamda#的退出

                                break#从反射点循环退出，控制从if self.condition_G(XR)的退出
                            else:
                                a0 *= lamda#反射点不可行
                    else:#中心点不可行
                        Xa = self.XL
                        Xb = Xcenter
                        for i in range(self.i):
                            self.X0_R[i][0] = Xa[i]
                            self.X0_R[i][1] = Xb[i]
                        self.Xk = self.random_X(self.k,[])
                        print(f'中心点不可行，更新随机区间={self.X0_R}')
                        break#控制else:#中心点不可行的退出
        print(f'迭代次数k={ci}')
xi = symbols('x1 x2')
F = (xi[0]-5)**2+4*(xi[1]-6)**2
g1 = 64-xi[0]**2-xi[1]**2
g2 = xi[1]-xi[0]-10
g3 = xi[0]-10
G = [g1,g2,g3]
X0_R = [[4,10],[5,15]]
k = 3
X0 = [7,9]
q1 = Complex_method(F,G,xi,X0_R,k,X0)
#收敛精度delta, 每次缩减的倍数lamda,反射a0 = 1.3, a0的最小限制a0_limit = e - 10，b0，扩张a1 = 1, b1, 缩小a2 = 0.5, b2, 压缩a3 = 0.5
q1.find_(1E-9,0.5,1,1E-10)