ActiveSetMethod.py

ActiveSetMethod.py

# -*- coding: utf-8 -*-
import numpy as np

def main():
    m = 2  # 制約条件の個数
    n = 2  # 設計変数の個数
    x0 = np.matrix(np.zeros([n, 1]))  # x初期値(適当に0ベクトルで与える)

    print("-" * 80 + "\n■2つある制約条件のうちどちらも無効制約条件の場合")
    Q = np.matrix([[2, 0], [0, 2]])
    c = np.matrix([[0], [0]])
    A = np.matrix([[1, 1], [-1, 1]])
    b = np.matrix([[1], [1]])
    opt = ActiveSetMethod(Q, c, A, b)
    opt.optLoop(x0)

    print("-" * 80 + "\n■2つある制約条件のうち1つが有効制約条件の場合(2番目の制約条件が有効制約)")
    Q = np.matrix([[2, 0], [0, 2]])
    c = np.matrix([[0], [0]])
    A = np.matrix([[1, 1], [-1, 1]])
    b = np.matrix([[1], [-1]])
    opt = ActiveSetMethod(Q, c, A, b)
    opt.optLoop(x0)

    print("-" * 80 + "\n■2つある制約条件のうちどちらも有効制約条件の場合")
    Q = np.matrix([[2, 0], [0, 2]])
    c = np.matrix([[0], [0]])
    A = np.matrix([[1, 1], [-1, 1]])
    b = np.matrix([[-1], [-1]])
    opt = ActiveSetMethod(Q, c, A, b)
    opt.optLoop(x0)

    print("-" * 80 + "\n■目的関数と制約関数を乱数で作成した場合")
    Q = np.matrix(np.random.randn(n, n))
    c = np.matrix(np.random.randn(n, 1))
    A = np.matrix(np.random.randn(m, n))
    b = np.matrix(np.random.randn(m, 1))
    opt = ActiveSetMethod(Q, c, A, b)
    opt.optLoop(x0)


class ActiveSetMethod(object):
    def __init__(self, Q, c, A, b, maxIter=1000):
        self.Q = Q
        self.c = c
        self.A = A
        self.b = b
        self.maxIter = maxIter  # 最適化ループの最大回数
        self.maxIterGetInitialFeasibleX = 1000  # 初期実行可能解探索ループの最大回数
        self.alpha = 1.1  # 初期実行可能解探索用の係数
        self.eps = 1e-10  # 0判別値

    def getG(self, x):
        """制約関数の値を取得"""
        return self.A * x - self.b

    def getF(self, x):
        """目的関数の値を取得"""
        return 0.5 * x.T * self.Q * x - self.c.T * x

    def getGradG(self):
        """制約関数の勾配を取得"""
        return self.A.T

    def getGradF(self, x):
        """目的関数の勾配を取得"""
        return self.Q * x - self.c

    def getSolution(self, activeSetIndexArray):
        """等式制約時の解とラグランジュ乗数を取得"""
        m, n = np.shape(self.A)
        if len(activeSetIndexArray) == 0:  # 制約条件が全て無効制約の場合
            xHat = np.linalg.inv(self.Q) * self.c
            lambdaHat = np.zeros(m)
        else:  # 有効制約条件がある場合
            Aq = self.A[activeSetIndexArray, :]
            bq = self.b[activeSetIndexArray]
            mq, nq = np.shape(Aq)
            tmp = np.linalg.inv(np.r_[np.c_[Aq, np.zeros([mq, mq])], np.c_[self.Q, Aq.T]]) * np.r_[bq, self.c]
            xHat = tmp[0:n]
            lambdaHat = np.zeros(m)
            lambdaHat[activeSetIndexArray] = tmp[n:len(tmp)]
        return xHat, np.matrix(lambdaHat).T

    def getActiveSetIndex(self, x):
        """有効制約条件番号を取得"""
        return np.where(np.abs(self.getG(x)) < self.eps)[0]

    def getInitialFeasibleX(self, x):
        """初期実行可能解の取得"""
        for iterGetInitialFeasibleX in range(self.maxIterGetInitialFeasibleX):
            maxGIndex = np.where(self.getG(x) == np.max(self.getG(x)))[0][0]
            x -= self.alpha * self.getGradG()[:, maxGIndex] * np.abs(self.A[maxGIndex] * x - self.b[maxGIndex]) / np.linalg.norm(self.A[maxGIndex]) ** 2
            if self.discriminationFeasiblity(x):
                break
        return x

    def discriminationFeasiblity(self, x):
        """実行可能解を判別"""
        return (sum(self.getG(x) < self.eps) == len(self.getG(x))).tolist()[0][0]

    def getT(self, xHat, x):
        """xHatを実行可能解に修正するときに使う係数tを取得"""
        gradGwrtT = self.A * (xHat - x)
        gradGwrtT[np.where(gradGwrtT == 0.0)] = self.eps
        T = - self.getG(x) / gradGwrtT
        if len(T[np.where(T > 0.0)].tolist()[0]) == 0:
            t = 0.0
        else:
            t = np.min(T[np.where(T > 0.0)])
        return t

    def getXHatDash(self, xHat, x):
        """xHatを実行可能解xHatDashに修正"""
        return x + self.getT(xHat, x) * (xHat - x)

    def getKKT(self, x, lambda_):
        """KKT条件を取得"""
        kkt1 = (self.getGradF(x) + self.getGradG() * lambda_).T.tolist()[0]
        kkt1Flag = sum(np.abs(kkt1) <= self.eps) == len(kkt1)
        kkt2 = self.getG(x).T.tolist()[0]
        kkt2Flag = sum(np.array(kkt2) <= self.eps) == len(kkt2)
        kkt3 = lambda_.T.tolist()[0]
        kkt3Flag = sum(np.array(kkt3) >= -self.eps) == len(kkt3)
        kkt4 = (np.diag(lambda_.T.tolist()[0]) * self.getG(x)).T.tolist()[0]
        kkt4Flag = sum(np.abs(kkt4) <= self.eps) == len(kkt4)
        return {"KKT1": [[kkt1Flag], kkt1], "KKT2": [[kkt2Flag], kkt2], "KKT3": [[kkt3Flag], kkt3],"KKT4": [[kkt4Flag], kkt4]}

    def step1(self, x0):
        """最適化ルーチンステップ1"""
        if self.discriminationFeasiblity(x0):
            x = x0
        else:
            x = self.getInitialFeasibleX(x0)
        return x, self.getActiveSetIndex(x)

    def step2(self, W):
        """最適化ルーチンステップ2"""
        return self.getSolution(W)

    def step3(self, xHat, x):
        """最適化ルーチンステップ3"""
        xHatDash = self.getXHatDash(xHat, x)
        return xHatDash, self.getActiveSetIndex(xHatDash)

    def step4(self, lambdaHat, W):
        """最適化ルーチンステップ4"""
        if sum(lambdaHat >= 0.0) == len(lambdaHat):
            optFlag = True
        else:
            minLambdaIndex = np.where(lambdaHat == np.min(lambdaHat))
            W = np.sort(np.delete(W, minLambdaIndex))
            optFlag = False
        return optFlag, W

    def optLoop(self, x0):
        """最適化ルーチンのループ"""
        x, W = self.step1(x0)

        if not self.discriminationFeasiblity(x):  # 初期実行可能解が見つからない場合
            print("初期実行可能解を見つけることができませんでした。x0の値を変更して下さい。")
            iter = 0
            optFlag = False
            return

        for iter in range(1, self.maxIter + 1):
            xHat, lambdaHat = self.step2(W)

            if not self.discriminationFeasiblity(xHat):
                x, W = self.step3(xHat, x)
                optFlag = False
            else:
                optFlag, W = self.step4(lambdaHat, W)

            self.xOpt_ = xHat
            self.lambdaOpt_ = lambdaHat

            if optFlag:  # 最適解が見つかった場合
                kkt = self.getKKT(self.xOpt_, self.lambdaOpt_)
                print(str(iter) + "回目のイタレーションで最適解が見つかりました。")
                print("・最適解x: " + str(self.xOpt_.T.tolist()[0]))
                print("・目的関数f(x)の最小値: " + str(self.getF(self.xOpt_).tolist()[0][0]))
                print("・KKT条件1(∇f+λ∇g=0): " + str(kkt["KKT1"][0][0]) + " " + str(kkt["KKT1"][1]))
                print("・KKT条件2(g<=0): " + str(kkt["KKT2"][0][0]) + " " + str(kkt["KKT2"][1]))
                print("・KKT条件3(λ>=0): " + str(kkt["KKT3"][0][0]) + " " + str(kkt["KKT3"][1]))
                print("・KKT条件4(λg=0): " + str(kkt["KKT4"][0][0]) + " " + str(kkt["KKT4"][1]))
                break

        if not optFlag:  # maxIter回数のループを回しても最適解が見つからなかった場合
            kkt = self.getKKT(self.xOpt_, self.lambdaOpt_)
            print(str(iter) + "回のイタレーションで最適解を見つけることができませんでした。")
            print("・最終イタレーションの設計変数x: " + str(self.xOpt_.T.tolist()[0]))
            print("・最終イタレーションの目的関数f(x)値: " + str(self.getF(self.xOpt_).tolist()[0][0]))
            print("・KKT条件1(∇f+λ∇g=0): " + str(kkt["KKT1"][0][0]) + " " + str(kkt["KKT1"][1]))
            print("・KKT条件2(g<=0): " + str(kkt["KKT2"][0][0]) + " " + str(kkt["KKT2"][1]))
            print("・KKT条件3(λ>=0): " + str(kkt["KKT3"][0][0]) + " " + str(kkt["KKT3"][1]))
            print("・KKT条件4(λg=0): " + str(kkt["KKT4"][0][0]) + " " + str(kkt["KKT4"][1]))

if __name__ == "__main__":
    main()