【运筹优化】粒子群算法(PSO)实现等误差直线逼近复杂曲线（附PSO Java源码）

本文链接：https://blog.youkuaiyun.com/weixin_52808620/article/details/134678537

本文介绍了传统几何等误差逼近算法及其在处理复杂曲线时的局限性，随后阐述了基于粒子群优化的等误差逼近算法，通过数学模型和粒子群求解，有效避免了传统方法中的问题。作者还提供了Java实现的粒子群算法示例。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

这篇文章是某工机电工程学院数控课设的产物，我实现了一个简单的逼近算法。

等误差算法要求每段逼近直线与自由曲线轮廓的误差相等，能够有效控制走刀路径的加工精度，减少逼近的直线段的数量。传统的等误差逼近采用几何关系计算节点坐标，这种方法适用于对曲线要求不高的场景，但如果曲线过于复杂，算法会面临编程复杂且容易出现失效的问题。为此，目前出现了一种使用基于粒子群优化算法进行等误差逼近的算法，可以适用于复杂曲线的逼近。下面对这两种等误差逼近算法进行介绍。

1、传统几何等误差逼近算法

使用传统几何等误差逼近的步骤如下：
（1）在起点A处作半径等于误差 $e$ 的圆方程，并将圆方程改写为 $x$ 关于 $y$ 的表达式 $(x)y=c\ \left( x \right)$ 。
（2）求误差圆与曲线的公切线PT的斜率 $k=yT−ypxT−xpk=\frac{y_T-y_p}{x_T-x_p}$ 。为了求得 $x_T,y_T,x_p,y_p$ ，需要联立以下方程组：
公差圆方程： $y=c(x)y=c\left( x \right)$
圆切线方程： $yT−yP=c′(XP)(xT−xP)y_T-y_P=c'\left( X_P \right) \left( x_T-x_P \right)$
曲线方程： $yT=f(xT)y_T=f\left( x_T \right)$
曲线在P点的切线方程： $yT−yP=f’(xT−xP)y_T-y_P=f’\left( x_T-x_P \right)$
（3）求弦长AB方程。使得AB平行于PT，即AB斜率为 $k$ ，得到AB的方程为： $y−ya=k(x−xa)y-y_a=k\left( x-x_a \right)$ 。
（4）联立曲线方程 $y=f(x)y=f\left( x \right)$ 与AB方程 $y−ya=k(x−xa)y-y_a=k\left( x-x_a \right)$ 求得B点的坐标。
（5）重复上述步骤求出曲线其余点坐标。
在这里插入图片描述

这种方法在曲线较复杂的情况下会出现失效。当非圆曲线出现局部形状变化异常或小变形时，误差过大可能会造成非圆曲线与直线段求交出现错误交点或无交点，如下图所示，曲线1出现了错误交点的情况，而曲线2出现了无交点的情况。
在这里插入图片描述

2、基于粒子群优化算法的等误差逼近算法

传统等误差逼近算法通过几何关系计算节点坐标，对曲线有一定要求。粒子群算法是一种群智能优化算法，将等误差逼近转化为一个优化数学模型，使用粒子群算法进行优化，可以有效解决传统几何等误差逼近中出现错误交点或者无交点的情况。

将自由曲线与逼近直线的误差作为目标函数，约束条件为逼近节点必须在自由曲线上，以一定的误差允许范围作为迭代终止条件构建出数学模型，通过粒子群算法求解数学模型即可得到下一个节点坐标。数学模型建立如下：

（1）设自由曲线方程为 $y=f(x)y=f\left( x \right)$ ，允许误差为 $e$ ，波动误差为 $b_n$ ，逼近直线段的起点坐标为 $(X0,Y0)\left( X_0,Y_0 \right)$ ，逼近直线段的终点坐标为 $(X,Y)\left( X,Y \right)$ ，自由曲线上的点的坐标为 $(xj,yj)\left( x_j,y_j \right)$ ，逼近直线的斜率为 $k=(Y−Y0)(X−X0)k=\frac{\left( Y-Y_0 \right)}{\left( X-X_0 \right)}$ ，逼近直线的方程为 $kx-y-kX_0+Y_0=0$ ，逼近直线与自由曲线之间的误差距离为 $d=\frac{|kx_j-y_j-kX_0+Y_0|}{\sqrt{1+k^2}}$
（2）目标函数为最小化误差距离与允许误差之间的差距，表示为： $F\left( x \right) =d-e$
（3）约束函数为所求的点必须在自由曲线上，表示为： $y=f\left( x \right)$
（4）迭代终止条件为： $F\left( x \right) \le b_n$
以上优化模型可以使用粒子群算法进行求解。粒子群优化算法(PSO：Particle swarm optimization) 是一种进化计算技术（evolutionary computation）。源于对鸟群捕食的行为研究。粒子群优化算法的基本思想：是通过群体中个体之间的协作和信息共享来寻找最优解。使用粒子群算法可以避免麻烦的几何计算和由于几何计算引起的算法失效。

基于粒子群算法优化后的等误差逼近的流程为：输入逼近误差、波动误差，配置粒子群算法的参数，之后在曲线逼近的范围中随机初始化各个函数值的粒子，评估每个粒子并得到全局最优。如果未至迭代结束条件，就不断更新粒子速度位置、计算每个粒子的函数值，让粒子自我学习以及社会学习，直到收敛至最优或者满足迭代结束条件，结束算法，输出逼近节点。由于算法不需要进行相切、相交等几何计算，所以不会出现上文所述的错误交点和无交点情况。
在这里插入图片描述

下面是我实现的利用pso求解连续函数的Java源码。稍微改进一下即可利用上述公式求出等误差逼近的节点。

实体类

package person.ghl.psoalgorithm.pojo;

import person.ghl.psoalgorithm.config.PsoConfig;

import java.util.Random;

/**
 * 粒子类
 * @author xcy
 */
public class Particle {
    /**
     * 粒子目前位置
     */
    private double[] nowPosition=new double[PsoConfig.dimension];
    /**
     * 粒子历史最优位置
     */
    private double[] bestPosition=new double[PsoConfig.dimension];
    /**
     * 粒子历史最优适应度
     */
    private double bestFitnessValue;
    /**
     * 粒子速度
     */
    private double[] velocity=new double[PsoConfig.dimension];
    /**
     * 粒子适应值
     */
    private double nowFitnessValue;
    /**
     * 随机数生成器
     */
    private final Random random=new Random();




    /**
     * 对速度、目前位置、最佳位置、最佳适应度进行初始化
     */
    public void initMyself() throws Exception {
        //速度随机在最大允许速度之内
        double doubleMaxVel=PsoConfig.maxVelocity*2;
        for(int i=0;i<this.velocity.length;i++){
            velocity[i]=random.nextDouble()*doubleMaxVel-PsoConfig.maxVelocity;
        }
        //目前位置要在最大位置允许范围之内
        for(int i=0;i<this.nowPosition.length;i++){
            nowPosition[i]=PsoConfig.maxPosition[i]-(PsoConfig.maxPosition[i]-PsoConfig.minPosition[i])*random.nextDouble();
            bestPosition[i]=nowPosition[i];
        }
        nowFitnessValue=PsoConfig.calculateFitnessValue(this);
        bestFitnessValue=nowFitnessValue;
    }

    /**
     * 更新粒子的当前位置
     * @param globalBestPosition:全局最优位置
     */
    public void updatePosAndFitness(double[] globalBestPosition, double velocityW) throws Exception {
        //下一代速度
        for(int i=0;i<this.velocity.length;i++){
            this.velocity[i]=velocityW*this.velocity[i]
                    + PsoConfig.c1*random.nextDouble()*(bestPosition[i]-nowPosition[i])
                    + PsoConfig.c2*random.nextDouble()*(globalBestPosition[i]-nowPosition[i]);
            //如果速度超过最大速度，限制最大速度
            if(Math.abs(velocity[i])>PsoConfig.maxVelocity){
                velocity[i]=velocity[i]<0?-PsoConfig.maxVelocity:PsoConfig.maxVelocity;
            }
        }
        //下一代位置
        for(int i=0;i<this.nowPosition.length;i++){
            nowPosition[i]=nowPosition[i]+velocity[i];
            //如果位置超过边界，限制在边界
            if(nowPosition[i]>PsoConfig.maxPosition[i]){
                nowPosition[i]=PsoConfig.maxPosition[i];
            }else if(nowPosition[i]<PsoConfig.minPosition[i]){
                nowPosition[i]=PsoConfig.minPosition[i];
            }
        }
        this.nowFitnessValue=PsoConfig.calculateFitnessValue(this);
    }


    public double[] getNowPosition() {
        return nowPosition;
    }

    public double[] getBestPosition() {
        return bestPosition;
    }

    public double getNowFitnessValue() {
        return nowFitnessValue;
    }

    public void setNowFitnessValue(double nowFitnessValue) {
        this.nowFitnessValue = nowFitnessValue;
    }

    public double getBestFitnessValue() {
        return bestFitnessValue;
    }

    public void setBestFitnessValue(double bestFitnessValue) {
        this.bestFitnessValue = bestFitnessValue;
    }

    public void setBestPosition(double[] bestPosition) {
        this.bestPosition = bestPosition;
    }

    public void setNowPosition(double[] nowPosition) {
        this.nowPosition = nowPosition;
    }

    public void setVelocity(double[] velocity) {
        this.velocity = velocity;
    }

    public Particle cloneMyself(){
        Particle clone=new Particle();
        clone.setBestPosition(this.bestPosition.clone());
        clone.setBestFitnessValue(this.bestFitnessValue);
        clone.setNowFitnessValue(this.nowFitnessValue);
        clone.setNowPosition(this.nowPosition.clone());
        clone.setVelocity(this.velocity.clone());
        return clone;
    }
}

配置类

package person.ghl.psoalgorithm.config;

import person.ghl.psoalgorithm.pojo.Particle;

import java.util.ArrayList;
import java.util.LinkedList;

/**
 * 粒子群算法的配置类
 * 如果是求最小值的话，只需要对目标函数以及变量数以及变量定义域进行修改就可以了
 * 如果是求最大值，除了修改上述内容以外，还需要对方法中的内容进行修改。
 *
 * @author xcy
 */
public class PsoConfig {
    /**
     * 速度权重的最大值和最小值，第一个是最小值，第二个是最大值
     */
    public static double[] w = {0.5, 0.9};
    /**
     * 自我学习因子
     */
    public static double c1 = 0.25;
    /**
     * 社会学习因子
     */
    public static double c2 = 0.35;
    /**
     * 最大迭代次数
     */
    public static int maxGeneration = 100;
    /**
     * 粒子群中粒子个数
     */
    public static int particleSize = 300;
    /**
     * 粒子位置上边界(需要限制每一个变量的范围，有多少个变量就在数组中写入多少个值)注意是上边界！！！不要写下边界在里面
     */
    public static double[] maxPosition;
    /**
     * 粒子位置下边界(需要限制每一个变量的范围)
     */
    public static double[] minPosition;
    /**
     * 粒子允许最大速度
     */
    public static double maxVelocity = 2.5;
    /**
     * 变量个数
     */
    public static int dimension = 1;

    /**
     * 一些常数/系数，根据函数形状自己增添系数
     */
    public static double[] beginPoint; //非圆曲线的起点的x坐标和y坐标
    /**
     * 允许误差
     */
    public static double e;
    /**
     * 波动误差
     */
    public static double bn;


    /**
     * 计算适应值（这个函数用来编辑满足约束条件的目标函数）
     */
    public static double calculateFitnessValue(Particle particle) throws Exception {
        //当前求出的解
        double[] variable = particle.getNowPosition();
        //算k
        double[] currentHeartLineValue = heartLine(variable[0]);
        double k = (currentHeartLineValue[1] - beginPoint[1]) / (currentHeartLineValue[0] - beginPoint[0]);
        //算适应度  从起始点到当前解搜索50次以特定步长找出最大
        double acceleration = (variable[0] - minPosition[0]) / 50;
        double bestLoss = Double.MIN_VALUE;
        double startT = minPosition[0];
        while (startT < variable[0]) {
            double[] heartLine = heartLine(startT);
            double currentLoss = calcLoss(k, heartLine[0], heartLine[1]);
            if (currentLoss > bestLoss) {
                bestLoss = currentLoss;
            }
            startT += acceleration;
        }
        return Math.abs(bestLoss - e);
    }

    /**
     * 计算误差
     */
    public static double calcLoss(double k, double xj, double yj) {
        return Math.abs(k * xj - yj - k * beginPoint[0] + beginPoint[1]) / Math.sqrt(1 + k * k);
    }

    /**
     * 心形线函数
     */
    public static double[] heartLine(double t) {
        double x = 40 * (2 * Math.cos(t) - Math.cos(2 * t));
        double y = 40 * (2 * Math.sin(t) - Math.sin(2 * t));
        return new double[]{x, y};
    }


    /**
     * 更新历史最优值(目标求最大值，if大于，目标求最小值，if小于)
     */
    public static void updateBestPosAndFitness(Particle particle) {
        //目前适应度值
        double nowFitnessValue = particle.getNowFitnessValue();
        //历史最优适应度值
        double bestFitnessValue = particle.getBestFitnessValue();
        if (nowFitnessValue < bestFitnessValue) {
            particle.setBestFitnessValue(nowFitnessValue);
            particle.setBestPosition(particle.getNowPosition().clone());
        }
    }

    /**
     * 返回最优粒子(需要根据具体情况选择升序还是降序，求最大值使用升序，求最小值使用降序)
     */
    public static Particle getBestParticle(ArrayList<Particle> particles) {
        LinkedList<Particle> particlesCloned = new LinkedList<>(particles);
        //对列表进行降序排序
        particlesCloned.sort((o1, o2) -> (int) (o2.getBestFitnessValue() * 1e20 - o1.getBestFitnessValue() * 1e20));
        return particlesCloned.getLast();
    }

    /**
     * 更新全局最优适应度以及全局最优位置(根据目标是最大，if是大于，目标是最小，if是小于)
     */
    public static double updateGlobalBestPosAndFitness(Particle bestParticle, double globalBestFitnessValue, double[] globalBestPosition) {
        if (bestParticle.getBestFitnessValue() < globalBestFitnessValue) {
            for (int i = 0; i < globalBestPosition.length; i++) {
                globalBestPosition[i] = bestParticle.getBestPosition()[i];
            }
            globalBestFitnessValue = bestParticle.getBestFitnessValue();
        }
        return globalBestFitnessValue;
    }

}

算法类

package person.ghl.psoalgorithm;



import person.ghl.psoalgorithm.config.PsoConfig;
import person.ghl.psoalgorithm.pojo.Particle;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.LinkedList;

/**
 * 粒子群算法实现类
 * @author xcy
 */
public class PsoAlgorithm {

    /**
     * 粒子群算法入口
     */
    public double solve() throws Exception {
        long begin=System.currentTimeMillis();
        //初始化粒子群
        ArrayList<Particle> particles=initParticles();
        //全局最优粒子位置以及其适应度
        double globalBestFitnessValue;
        double[] globalBestPosition;
        Particle bestParticle= PsoConfig.getBestParticle(particles);
        globalBestPosition=bestParticle.getBestPosition().clone();
        globalBestFitnessValue=bestParticle.getBestFitnessValue();

        for(int i=0;i<PsoConfig.maxGeneration;i++){
            //更新粒子位置和适应值，所有粒子进入下一代
            //速度权重
            double velocityW= ((PsoConfig.w[1]- PsoConfig.w[0])
                    *(PsoConfig.maxGeneration-i)/ PsoConfig.maxGeneration)
                    + PsoConfig.w[0];
            for(Particle particle:particles){
                particle.updatePosAndFitness(globalBestPosition, velocityW);
            }
            //更新历史最优值
            for(Particle particle:particles){
                PsoConfig.updateBestPosAndFitness(particle);
            }
            //更新全局最优值和全局最优位置
            bestParticle=PsoConfig.getBestParticle(particles).cloneMyself();
            globalBestFitnessValue=PsoConfig.updateGlobalBestPosAndFitness(bestParticle,globalBestFitnessValue,globalBestPosition);
            //这里根据需要插入一个迭代终止条件
            if(globalBestFitnessValue<=PsoConfig.bn){
                break;
            }
        }
        //输出最优位置和最优值
        System.out.println("------------------------");
        System.out.println("最优位置："+ Arrays.toString(globalBestPosition));
        System.out.println("最优适应度（误差-允许误差）："+globalBestFitnessValue);
        System.out.println("误差："+(globalBestFitnessValue+PsoConfig.e));
        System.out.println("求解时间："+(System.currentTimeMillis()-begin)+" ms");
        System.out.println("------------------------");
        return globalBestPosition[0];
    }

    /**
     * 初始化粒子群
     */
    private ArrayList<Particle> initParticles() throws Exception {
        LinkedList<Particle> particles= new LinkedList<>();
        for(int i = 0; i< PsoConfig.particleSize; i++){
            Particle particle=new Particle();
            particle.initMyself();
            particles.add(particle);
        }
        return new ArrayList<>(particles);
    }






}