使用粒子群算法求解约束优化问题

一、使用罚函数法将约束优化问题转化为无约束优化问题

无约束优化问题，是众多优化问题中最基本的问题，它对自变量$x$的取值范围不加限制，多以无需考虑$x$的可行性。
$$\begin{array}{r}\underset{x\in {\mathbb{R}}^n}{\min }f(x)\end{array}$$

约束优化问题中，自变量$x$不能任意取值这导致许多无约束优化算法不能使用。 例如梯度法中沿着负梯度方向下降所得的点未必是可行点，要寻找的最优解处目标函数的梯度也不是零向量，这使得约束优化问题比无约束优化问题要复杂许多。
$$\begin{array}{rl}\min & f(x)\\ \text{s.t}& x\in \mathcal{X}\end{array}$$

罚函数法（Penalty Function Approach）
求解约束优化问题最常采用的方法是罚函数法（Penalty Function Approach），将约束作为惩罚项加到目标函数中，从而转化为我们熟悉的无约束优化问题求解。

等式约束的二次罚函数

$$\begin{array}{rl}\min & f(x)\\ \text{s.t}& c_i(x)=0,i\in \epsilon \end{array}$$
其中,$x\in {\mathbb{R}}^n$，$\epsilon$为等式约束的指标集，$c_i(x)$为连续函数．

罚函数法的思想是将约束优化问题转化为无约束优化问题来进行求解。为了保证解的逼近质量，无约束优化问题的目标函数为原约束优化问题的目标函数加上与约束函数有关的惩罚项。对于可行域外的点，惩罚项为正，即对该点进行惩罚；对于可行域内的点，惩罚项为0，即不做任何惩罚。因此，惩罚项会促使无约束优化问题的解落在可行域内。对于等式约束问题，惩罚项的选取方式有很多，结构最简单的是二次函数。这里给出二次罚函数的定义：
$$\begin{array}{r}{P}_E(x,\sigma )=f(x)+\frac{1}{2}\sigma \sum _{i\in \epsilon }{c}_i^2(x)\end{array}$$
其中等式右端第二项称为惩罚项，$\sigma >0$称为罚因子。

由于这种罚函数对不满足约束的点进行惩罚，在迭代过程中点列一般 处于可行域之外，因此它也被称为外点罚函数。二次罚函数的特点如下：对于非可行点而言，当σ变大时，惩罚项在罚函数中的权重加大，对罚函数求极小，相当于迫使其极小点向可行域靠近；在可行域中，$P_E(x,\sigma )$的全局极小点与约束最优化问题(4-5)的最优解相同。

$$\begin{array}{rl}\min & x+\sqrt{3}y\\ \text{s.t}& x^2+y^2=1\end{array}$$
该问题的最优解为$(-\frac{1}{2},-\frac{\sqrt{3}}{2})$。

考虑二次罚函数：
$$P_E(x,y,\sigma )=x+\sqrt{3}y+\frac{1}{2}\sigma (x^2+y^2-1{)}^2$$

不等式约束优化问题的二次罚函数法

$$\begin{array}{rl}\min & f(x)\\ \text{s.t}& c_i(x)\le 0,i\in \epsilon \end{array}$$
显然，它和等式约束优化问题最大的不同就是允许 $c_i(x)<0$ 发生，而若采用原来的方式定义罚函数为 $\Vert c(x){\Vert }^2$，它也会惩罚 $c_i(x)<0$ 的可行点，这显然不是我们需要的，必须对原有二次罚函数进行改造来得到新的二次罚函数，它应该具有如下特点：仅仅惩罚$c_i(x)>0$的那些点，而对可行点不作惩罚：

$$P_I(x,\sigma )=f(x)+\frac{1}{2}\sigma \sum _{i\in I}{\tilde{c}}_i^2(x)$$
其中等式右端第二项称为惩罚项， c ~ i 2 ( x ) = max ⁡ { c i ( x ) , 0 } \tilde{c}_i^2(x)=\max\{c_i(x), 0\} c~i2​(x)=max{ci​(x),0}。

注意到函数 h ( t ) = ( max ⁡ { t , 0 } ) 2 h(t) = (\max\{t, 0\})^2 h(t)=(max{t,0})2是关于$t$可导的，因此$P_I(x,\sigma )$的梯度也存在，可以使用梯度类算法来求解子问题。

一般约束问题的二次罚函数法

一般的约束优化问题可能既含等式约束又含不等式约束，它的形式为：
$$\begin{array}{rl}\min & f(x)\\ \text{s.t}& c_i(x)=0,i\in \epsilon \\ & c_i(x)\le 0,i\in \mathcal{I}\end{array}$$
针对这个问题，我们只需要将两种约束的罚函数相加就能得到一般约束优化问题的二次罚函数：

$$P(x,\sigma )=f(x)+\frac{1}{2}\sigma [\sum _{i\in \epsilon }{c}_i^2(x)+\sum _{i\in I}{\tilde{c}}_i^2(x)]$$
其中:

• 等式右端第二项称为惩罚项， c ~ i 2 ( x ) = max ⁡ { c i ( x ) , 0 } \tilde{c}_i^2(x)=\max\{c_i(x), 0\} c~i2​(x)=max{ci​(x),0}。
• $\sigma$为惩罚因子

对参数${\sigma }_k$的选取需要非常小心，若${\sigma }_k$增长太快，则子问题不易求解（具体分析见下一小节末尾对算法数值困难的讨论）。若增长太慢，则算法需要的外迭代数（算法中的 while 循环）会增多。一个比较合理的取法是根据当前$P_E(x,{\sigma }_k)$的求解难度来确定${\sigma }_k$的增幅，若当前子问题收敛很快，可以在下一步选取较大的${\sigma }_{k+1}$，否则就不宜过分增大${\sigma }_k$．

二、基于粒子群算法的约束优化问题求解

算例1

$$\begin{array}{rl}\min & x+\sqrt{3}y\\ \text{s.t}& x^2+y^2=1\end{array}$$
该问题的最优解为$(-\frac{1}{2},-\frac{\sqrt{3}}{2})$。

考虑二次罚函数：
$$P_E(x,y,\sigma )=x+\sqrt{3}y+\frac{1}{2}\sigma (x^2+y^2-1{)}^2$$

import matplotlib.pyplot as plt
import numpy as np


class ParticleSwarmOptimization:
    def __init__(self):
        self.w = 1  # 惯性权重
        self.c1 = 2  # 加速系数
        self.c2 = 2  # 加速系数
        self.dimension = 2
        self.pop_size = 100
        self.max_iteration = 1000
        self.sigma = 10
        self.max_velocity = .2

    def evaluate(self, X):
        """
        评价函数: 二次罚函数
        :return: float
        """
        if X.ndim == 1:
            x = X[0]
            y = X[1]
        else:
            x = X[:, 0]
            y = X[:, 1]
        return x + np.sqrt(3) * y + 0.5 * self.sigma * (x ** 2 + y ** 2 - 1) ** 2

    def update_velocity(self, X, V, pbest, gbest):
        r1 = np.random.random(size=(self.pop_size, self.dimension))
        r2 = np.random.random(size=(self.pop_size, self.dimension))
        V = self.w * V + self.c1 * r1 * (pbest - X) + self.c2 * r2 * (gbest - X)

        V[V < -self.max_velocity] = np.random.uniform(low=-.2, high=.2)
        V[V > self.max_velocity] = np.random.uniform(low=-.2, high=.2)

        return V

    def update_position(self, X, V):
        X = X + V
        return X

    def run(self):
        trace = []
        # 初始化粒子的位置和速度
        X = np.random.uniform(low=-1, high=1, size=(self.pop_size, self.dimension))
        V = np.random.uniform(low=-0.1, high=0.1, size=(self.pop_size, self.dimension))
        X_fitness = self.evaluate(X)
        pbest = X
        gbest = X[np.argmin(X_fitness)]
        pbest_fitness = X_fitness
        gbest_fitness = self.evaluate(gbest)

        for k in range(self.max_iteration):
            V = self.update_velocity(X, V, pbest, gbest)  # 更新速度
            X = self.update_position(X, V)  # 更新位置
            X_fitness = self.evaluate(X)
            bst_fitness = np.min(X_fitness)

            # 更新粒子历史最优位置
            for j in range(self.pop_size):
                if X_fitness[j] < pbest_fitness[j]:
                    pbest[j] = X[j]
                    pbest_fitness[j] = X_fitness[j]
            # 更新群体的最优位置
            if bst_fitness < gbest_fitness:
                gbest = X[np.argmin(X_fitness)]
                gbest_fitness = bst_fitness
            trace.append((gbest, gbest_fitness))
        f = lambda x, y: x + np.sqrt(3) * y  # 目标函数
        print(gbest, f(*gbest))
        self.plot_objective_value(trace)

    def plot_objective_value(self, trace):
        obj_list = [obj for gbest, obj in trace]
        plt.plot(np.arange(self.max_iteration), obj_list, color="#191970", linewidth=2., alpha=1.)
        plt.grid(True)
        plt.show()


if __name__ == '__main__':
    # 最优解（-0.5, -0.8660254037844386）, 目标函数值为-2
    solution = np.array([-0.5, -0.8660254037844386])
    pso = ParticleSwarmOptimization()
    pso.run()
    # print(pso.objective_func([-0.5, -0.8660254037844386]))

算例2

$$\begin{array}{rl}\min & f(x)=(x_1-2{)}^2+(x_2-1{)}^2\\ s.t.& x_1=2x_2-1\\ & \frac{x_1^2}{4}+x_2^2-1\le 0\end{array}$$
该问题最优解为1.3934651。

一般约束的二次罚函数为：

min ⁡ f ( x ) = ( x 1 − 2 ) 2 + ( x 2 − 1 ) 2 + 1 2 σ [ ( x 1 − 2 x 2 + 1 ) 2 + ( max ⁡ { x 1 2 4 + x 2 2 − 1 , 0 } ) 2 ] \begin{align} \min \quad & f(x)=(x_1-2)^2+(x_2-1)^2 +\frac{1}{2}\sigma\left[ (x_1-2x_2+1)^2 + (\max\{\frac{x_1^2}{4}+x_2^2-1, 0\})^2\right]\\ \end{align} min​f(x)=(x1​−2)2+(x2​−1)2+21​σ[(x1​−2x2​+1)2+(max{4x12​​+x22​−1,0})2]​​

import matplotlib.pyplot as plt
import numpy as np


class Problem:
    def __init__(self):
        pass

    def func(self, variables) -> float:
        pass


class ParticleSwarmOptimization:
    def __init__(self):
        self.w = 1  # 惯性权重
        self.c1 = 2  # 加速系数
        self.c2 = 2  # 加速系数
        self.dimension = 2
        self.pop_size = 100
        self.max_iteration = 1000
        self.sigma = 100
        self.max_velocity = .2

    def evaluate(self, X):
        """
        评价函数: 二次罚函数
        :return: float
        """
        if X.ndim == 1:
            x = X[0]
            y = X[1]
            func2 = (x - 2) ** 2 + (y - 1) ** 2 + 0.5 * self.sigma * ((x - 2 * y + 1) ** 2 + max(x ** 2 / 4 + y ** 2 - 1, 0))
        else:
            x = X[:, 0]
            y = X[:, 1]
            columns = x ** 2 / 4 + y ** 2 - 1
            for i in range(len(columns)):
                columns[i] = max(0, columns[i])

            func2 = (x - 2) ** 2 + (y - 1) ** 2 + 0.5 * self.sigma * ((x - 2 * y + 1) ** 2 + columns)
        # func1 = x + np.sqrt(3) * y + 0.5 * self.sigma * (x ** 2 + y ** 2 - 1) ** 2
        return func2

    def update_velocity(self, X, V, pbest, gbest):
        r1 = np.random.random(size=(self.pop_size, self.dimension))
        r2 = np.random.random(size=(self.pop_size, self.dimension))
        V = self.w * V + self.c1 * r1 * (pbest - X) + self.c2 * r2 * (gbest - X)

        V[V < -self.max_velocity] = np.random.uniform(low=-.2, high=.2)
        V[V > self.max_velocity] = np.random.uniform(low=-.2, high=.2)

        return V

    def update_position(self, X, V):
        X = X + V
        return X

    def run(self):
        trace = []
        # 初始化粒子的位置和速度
        X = np.random.uniform(low=-1, high=1, size=(self.pop_size, self.dimension))
        V = np.random.uniform(low=-0.1, high=0.1, size=(self.pop_size, self.dimension))
        X_fitness = self.evaluate(X)
        pbest = X
        gbest = X[np.argmin(X_fitness)]
        pbest_fitness = X_fitness
        gbest_fitness = self.evaluate(gbest)

        for k in range(self.max_iteration):
            V = self.update_velocity(X, V, pbest, gbest)  # 更新速度
            X = self.update_position(X, V)  # 更新位置
            X_fitness = self.evaluate(X)
            bst_fitness = np.min(X_fitness)

            # 更新粒子历史最优位置
            for j in range(self.pop_size):
                if X_fitness[j] < pbest_fitness[j]:
                    pbest[j] = X[j]
                    pbest_fitness[j] = X_fitness[j]
            # 更新群体的最优位置
            if bst_fitness < gbest_fitness:
                gbest = X[np.argmin(X_fitness)]
                gbest_fitness = bst_fitness
            trace.append((gbest, gbest_fitness))
        f = lambda x, y: (x - 2) ** 2 + (y - 1) ** 2  # 目标函数
        print(gbest, f(*gbest))
        self.plot_objective_value(trace)

    def plot_objective_value(self, trace):
        obj_list = [obj for gbest, obj in trace]
        plt.plot(np.arange(self.max_iteration), obj_list, color="#191970", linewidth=2., alpha=1.)
        plt.grid(True)
        plt.show()


if __name__ == '__main__':
    # 最优解（-0.5, -0.8660254037844386）, 目标函数值为-2
    pso = ParticleSwarmOptimization()
    pso.run()

参考：

• 粒子群算法求解带约束优化问题 - 知乎 (zhihu.com)
• 粒子群算法求解带约束优化问题 源码实现_带约束条件的 粒子群-CSDN博客
• Parsopoulos K E, Vrahatis M N. Particle Swarm Optimization Method for Constrained Optimization Problems[J]. 2002.
• 刘浩洋、户将.最优化计算方法