最优化交替方向乘子算法python实现

from sympy import *
import random
import sys


class AltDir:
    def __init__(self, x, y, f1, f2, A1, A2, b, k, m):
        # x 是x组成的字符串,构造函数将其处理为数学符号组成的列向量
        # y 是y组成的字符串
        # f1 是f1(x)的字符串表达式，f2是f(y)组成的字符串表达式
        # A1、A2是列表，b是b组成的列表，构造函数会将其转化为列向量
        self.x = x
        self.y = y
        self.Trans_Str_To_Var()
        # 将其转化为列向量
        self.x = Matrix(self.x)
        self.y = Matrix(self.y)
        # 下面这几个句子为什么写成函数来调用不行呢
        for i in range(1, self.x.shape[0] + 1):
            # 次句甚妙
            exec('x{0} = self.x[i-1, 0]'.format(i))
        for i in range(1, self.y.shape[0] + 1):
            exec("y{0} = self.y[i-1, 0]".format(i))
        self.f1 = eval(f1)
        self.f2 = eval(f2)
        self.A1 = Matrix(A1)
        self.A2 = Matrix(A2)
        self.b = Matrix(b)
        # k表示拉格朗日乘子
        self.k = Matrix(k)
        self.m = m
        self.boj_fun = self.f1 + self.f2
        self.LagAuFun = 0
        # self.Creat_LagAuFun()
        # 对self.x和self.y进行初值赋值
        self.x_x = self.x
        self.y_y = self.y
        self.x_x = ones(self.x.shape[0], 1)
        self.y_y = ones(self.y.shape[0], 1)
        # 存储偏导数，两种情况的x,y的偏导
        self.par = [0, 0]
        # 记录目标函数的值
        self.Goal = []
        # 存储拉格朗日增广函数的赋值计算结果
        self.ValueLagAuFun = []
        # 存储每次x和y的变化情况
        self.x_change = []
        self.y_change = []

    # 将向量进行转置，该方法需要传入一个列表
    @staticmethod
    def Trans_List_To_Cvector(r_list):
        np_r_list = Matrix(r_list)
        c_vec = np_r_list.reshape(-1, 1)
        return c_vec

    # 将字符串转化为数学变量
    def Trans_Str_To_Var(self):
        var_x = self.x.split(' ')
        var_y = self.y.split(' ')
        for i in range(0, len(var_x)):
            var_x[i] = symbols(var_x[i], real=True)
        for i in range(0, len(var_y)):
            var_y[i] = symbols(var_y[i], real=True)
        self.x = var_x
        self.y = var_y

    # 构建目标函数的拉格朗日增广函数
    def Creat_LagAuFun(self):
        #        LagAuFun1 = self.f1 + self.f2 + self.k.reshape(1, -1) * (
        #                self.A1 * self.x + self.A2 * self.y - self.b) + m / 2 * (self.A1 * self.x + self.A2 * self.y - self.b)
        #        LagAuFun2 = self.f1 + self.f2 + self.k.reshape(1, -1) * (
        #                self.A1 * self.x + self.A2 * self.y - self.b) - m / 2 * (self.A1 * self.x + self.A2 * self.y - self.b)
        #        self.LagAuFun[0] = LagAuFun1
        #        self.LagAuFun[1] = LagAuFun2
        LagAuFun = self.f1 + self.f2 + (self.k.transpose() * (
                self.A1 * self.x + self.A2 * self.y - self.b))[0, 0] + m / 2 * (
                       (self.A1 * self.x + self.A2 * self.y - self.b).norm()) ** 2
        self.LagAuFun = LagAuFun

    # 对函数进行求偏导
    def GetPartical(self):
        par_x = diff(self.LagAuFun, self.x)
        par_y = diff(self.LagAuFun, self.y)
        self.par[0] = par_x
        self.par[1] = par_y

    # 对拉格朗日增广函数中的变量进行赋值计算
    def CaculateLagAuFun(self):
        for i in range(1, self.x.shape[0] + 1):
            exec('x{0} = self.x_x[i-1, 0]'.format(i))
        for i in range(1, self.y.shape[0] + 1):
            exec("y{0} = self.y_y[i-1, 0]".format(i))
        repx = [(self.x[i, 0], self.x_x[i, 0]) for i in range(0, self.x.shape[0])]
        repy = [(self.y[i, 0], self.y_y[i, 0]) for i in range(0, self.y.shape[0])]
        f1 = eval(str(self.f1))
        f2 = eval(str(self.f2))
        repf = [(self.f1, f1), (self.f2, f2)]
        rep = repx + repy + repf
        ValueLagAuFun = self.LagAuFun.subs(rep)
        self.ValueLagAuFun.append(ValueLagAuFun)

    # 对除偏导以外的变量进行赋值
    def AssignmentX(self, k=None):
        for i in range(1, self.x.shape[0] + 1):
            exec('x{0} = self.x_x[i-1, 0]'.format(i))
        for i in range(1, self.y.shape[0] + 1):
            exec("y{0} = self.y_y[i-1, 0]".format(i))
        repx = [(self.x[i, 0], self.x_x[i, 0]) for i in range(0, self.x.shape[0]) if i != k]
        repy = [(self.y[i, 0], self.y_y[i, 0]) for i in range(0, self.y.shape[0])]
        exec('x{} = symbols("x{}", real=True)'.format(k + 1, k + 1))
        f1 = eval(str(self.f1))
        exec('x{} = self.x_x[k, 0]'.format(k + 1))
        f2 = eval(str(self.f2))
        repf = [(self.f1, f1), (self.f2, f2)]
        rep = repx + repy + repf
        return rep

    def AssignmentY(self, k=None):
        for i in range(1, self.x.shape[0] + 1):
            exec('x{0} = self.x_x[i-1, 0]'.format(i))
        for i in range(1, self.y.shape[0] + 1):
            exec("y{0} = self.y_y[i-1, 0]".format(i))
        repx = [(self.x[i, 0], self.x_x[i, 0]) for i in range(0, self.x.shape[0])]
        repy = [(self.y[i, 0], self.y_y[i, 0]) for i in range(0, self.y.shape[0]) if i != k]
        f1 = eval(str(self.f1))
        exec('y{} = symbols("y{}", real=True)'.format(k + 1, k + 1))
        f2 = eval(str(self.f2))
        exec('y{} = self.y_y[k, 0]'.format(k + 1))
        repf = [(self.f1, f1), (self.f2, f2)]
        rep = repx + repy + repf
        return rep

    # 对x、y、k进行更新
    def Update(self):
        REPX = []
        REPY = []
        # 先对x变量进行处理
        for i in range(0, self.x.shape[0]):
            repx = self.AssignmentX(i)
            REPX.append(repx)

        for i in range(0, self.y.shape[0]):
            repy = self.AssignmentY(i)
            REPY.append(repy)

        ParValueX = []
        ParValueY = []
        for i in range(0, len(REPX)):
            value = solve(diff(self.LagAuFun.subs(REPX[i])))
            ParValueX.append(value)

        for i in range(0, len(REPY)):
            value = solve(diff(self.LagAuFun.subs(REPY[i])))
            ParValueY.append(value)
        print()
        # x的更新
        MinParX = []
        for i in range(0, len(ParValueX)):
            if len(ParValueX[i]) > 1:
                Val = []
                for j in range(0, len(ParValueX[i])):
                    v = self.LagAuFun.subs(REPX[i])
                    exec('x{} = ParValueX[i][j]'.format(i))
                    v = v.subs(eval('x{}'.format(i)), ParValueX[i][j])
                    Val.append(v)
                I_ndex = Val.index(min(Val))
                MinParX.append(ParValueX[i][I_ndex])
            else:
                if len(ParValueX[i]) == 0:
                    print('input error')
                    sys.exit()
                else:
                    MinParX.append(ParValueX[i][0])
                # MinParX.append(ParValueX[i][0])
        x_change = MinParX
        self.x_change.append(x_change)
        self.x_x = Matrix(MinParX)

        # y的更新

        MinParY = []
        for i in range(0, len(ParValueY)):
            if len(ParValueY[i]) > 1:
                Val = []
                for j in range(0, len(ParValueY[i])):
                    v = self.LagAuFun.subs(REPY[i])
                    exec('y{} = ParValueY[i][j]'.format(i))
                    v = v.subs(eval('y{}'.format(i)), ParValueY[i][j])
                    Val.append(v)
                I_ndex = Val.index(min(Val))
                MinParY.append(ParValueY[i][I_ndex])
            else:
                if len(ParValueY[i]) == 0:
                    print('input error')
                    sys.exit()
                else:
                    MinParY.append(ParValueY[i][0])
                # MinParY.append(ParValueY[i][0])

        y_change = MinParY
        self.y_change.append(y_change)
        self.y_y = Matrix(MinParY)

        # 更新k
        t = random.uniform(0.1, (1 + 5 ** 0.5) / 2 - 0.1)
        self.k = self.k + t * self.m * (self.A1 * self.x_x + self.A2 * self.y_y - self.b)
        print()

    # 进行迭代
    def Iter(self):
        # 构建拉格朗日增广函数
        self.Creat_LagAuFun()
        # 赋值进行计算
        self.CaculateLagAuFun()
        # 求偏导
        self.GetPartical()
        # 更新
        self.Update()
        k = 1
        while True:
            self.Creat_LagAuFun()
            self.GetPartical()
            self.CaculateLagAuFun()
            if Abs(self.ValueLagAuFun[k] - self.ValueLagAuFun[k - 1]) < 0.0001:
                break
            else:
                self.Update()
                k += 1
            print(self.ValueLagAuFun[-1])
        return self.ValueLagAuFun


# 使用eval来获取字符串表达式结果
# ------------------------------------------------------------------------------------
# x = 'x1 x2'  # n维度
# y = 'y1 y2'  # m维度
# f1 = 'x1 ** 2 + x2 ** 2'
# f2 = 'y1 ** 2 + y2 ** 2'
# A1 = [[1, 2]]  # pxn
# A2 = [[1, 2]]  # pxm
# b = [5]  # p
# k = [2]  # p
# m = 0.01
# --------------------------------------------------------------------------------------

# x = 'x1 x2'  # n维度
# y = 'y1 y2'  # m维度
# f1 = 'x1 ** 2 + x2 ** 2'
# f2 = 'y1 ** 2 + y2 ** 2'
# A1 = [[1, 2], [7, 6]]  # pxn
# A2 = [[1, 2], [6, 3]]  # pxm
# b = [5, 12]  # p
# k = [2, 3]  # p
# m = 0.01

# --------------------------------------------------------------------------------------

# x = 'x1 x2 x3'  # n维度
# y = 'y1 y2 y3'  # m维度
# f1 = 'x1 ** 2 + x2 ** 2 + x3'


# f2 = 'y1 ** 2 + y2 ** 2 + y3**2'
# A1 = [[1, 2, 5], [7, 6, 9]]  # pxn
# A2 = [[1, 2, 8], [6, 3, 15]]  # pxm
# b = [5, 12]  # p
# k = [2, 3]  # p
# m = 0.0001


#---------------------------------------------------------------------------------------

# x = 'x1 x2'  # n维度
# y = 'y1 y2 y3'  # m维度
# f1 = 'x1 ** 2 + x2 ** 2'
# f2 = 'y1 ** 2 + y2 ** 2 + y3**2'
# A1 = [[1, 2], [7, 6]]  # pxn
# A2 = [[1, 2, 8], [6, 3, 15]]  # pxm
# b = [5, 12]  # p
# k = [2, 3]  # p
# m = 0.001

#----------------------------------------------------------------------------------------

# x = 'x1 x2 x3'  # n维度
# y = 'y1 y2'  # m维度
# f1 = 'x1 ** 2 + x2 ** 2 + x3'
# f2 = 'y1 ** 2 + y2 ** 2'
# A1 = [[1, 2, 5], [7, 6, 9]]  # pxn
# A2 = [[1, 2], [6, 3]]  # pxm
# b = [5, 12]  # p
# k = [2, 3]  # p
# m = 0.0001

# ----------------------------------------------------------------------------------------
x = 'x1'
y = 'y1'
f1 = 'x1**2'
f2 = 'y1**2'
A1 = [1]
A2 = [1]
b = [5]
k = [2]
m = 2
# ---------------------------------------------------------------------------------------
AD = AltDir(x, y, f1, f2, A1, A2, b, k, m)
result = AD.Iter()
print(result)