数学建模常用算法:粒子群算法(PSO)求解二元函数最小值+限定x,y范围测试【java实现--详细注释+Matlab绘制粒子群飞行过程】

代码

package com.dam.heuristic.pso.test;

import java.util.List;
import java.util.Random;

public class PsoApi {

    //粒子数量
    private int particleNum;
    //个体学习因子,设置得越大,粒子越容易根据自己的想法飞行,若设置过大,容易跳出局部最优,但收敛较慢
    private double c1;
    //社会学习因子,设置得越大,粒子越容易根据群体的想法飞行,若设置过大,容易陷入局部最优,收敛较快
    private double c2;
    //速度最大值
    private double vMax;
    //速度的惯性权重
    private double w;
    //迭代次数
    private int genMax;

    public PsoApi(int particleNum, double c1, double c2, double vMax, double w, int genMax) {
        this.particleNum = particleNum;
        this.c1 = c1;
        this.c2 = c2;
        this.vMax = vMax;
        this.w = w;
        this.genMax = genMax;
    }

    /**
     * 求解
     */
    public double[][][] solve() {
        变量声明
        //存储粒子
        Particle[] particleArr;
        //所有粒子找到的最优解(由于问题为最小化问题,设置初始最优值为较大的数)
        double gBest = Double.MAX_VALUE;
        //群体最优解对应的x和y
        double bestX = 0;
        double bestY = 0;
        //随机数工具
        Random random = new Random();
        long start = System.currentTimeMillis();
        //存储每一代粒子所在位置
        double[][][] positionArr = new double[this.genMax][this.particleNum][2];

        初始化粒子
        particleArr = new Particle[this.particleNum];
        for (int i = 0; i < particleArr.length; i++) {
            //初始化粒子群,注意:这里设置每个粒子的速度一样,读者可以根据自己的喜爱进行设置
            //初始化粒子的坐标和速度
            particleArr[i] = new Particle(-1000, 1000, -1000, 1000, 0.01, 0.01, random);
            //初始化粒子的函数值(粒子还没有开始飞,当前位置肯定是找到过的最优位置啦)
            particleArr[i].setBestX(particleArr[i].getX());
            particleArr[i].setBestY(particleArr[i].getY());
            double pValue = this.objectFunction(particleArr[i].getX(), particleArr[i].getY());
            particleArr[i].setpBest(pValue);
            //由于问题为最小化问题,目标函数越小越好
            if (pValue < gBest) {
                bestX = particleArr[i].getX();
                bestY = particleArr[i].getY();
                gBest = pValue;
            }
        }

        开始求解
        for (int i = 0; i < this.genMax; i++) {
            //对每个粒子进行操作
            for (int j = 0; j < this.particleNum; j++) {
                ///更新速度
                //更新x轴方向上的速度
                particleArr[j].setxV(this.w * particleArr[j].getxV()
                        + this.c1 * random.nextDouble() * (particleArr[j].getBestX() - particleArr[j].getX())
                        + this.c2 * random.nextDouble() * (bestX - particleArr[j].getX()));
                //处理越界
                if (particleArr[j].getxV() > this.vMax) {
                    particleArr[j].setxV(this.vMax);
                }else if (particleArr[j].getxV() < -this.vMax) {
                    particleArr[j].setxV(-this.vMax);
                }

                //更新y轴方向上的速度
                particleArr[j].setyV(this.w * particleArr[j].getyV()
                        + this.c1 * random.nextDouble() * (particleArr[j].getBestY() - particleArr[j].getY())
                        + this.c2 * random.nextDouble() * (bestY - particleArr[j].getY()));
                //处理越界
                if (particleArr[j].getyV() > this.vMax) {
                    particleArr[j].setyV(this.vMax);
                } else if (particleArr[j].getyV() < -this.vMax) {
                    particleArr[j].setyV(-this.vMax);
                }

                ///更新位置
                double nextX = particleArr[j].getX() + particleArr[j].getxV();
                //处理越界
                if (nextX > particleArr[j].getxMax()) {
                    nextX = particleArr[j].getxMax();
                } else if (nextX < particleArr[j].getxMin()) {
                    nextX = particleArr[j].getxMin();
                }
                particleArr[j].setX(nextX);
                double nextY = particleArr[j].getY() + particleArr[j].getyV();
                //处理越界
                if (nextY > particleArr[j].getyMax()) {
                    nextY = particleArr[j].getyMax();
                } else if (nextY < particleArr[j].getyMin()) {
                    nextY = particleArr[j].getyMin();
                }
                particleArr[j].setY(nextY);

                ///更新粒子历史最优解和粒子全体最优解
                double pValue = this.objectFunction(particleArr[j].getX(), particleArr[j].getY());
                if (pValue < particleArr[j].getpBest()) {
                    particleArr[j].setBestX(particleArr[j].getX());
                    particleArr[j].setBestY(particleArr[j].getY());
                    particleArr[j].setpBest(pValue);
                }
                //由于问题为最小化问题,目标函数越小越好
                if (pValue < gBest) {
                    bestX = particleArr[j].getX();
                    bestY = particleArr[j].getY();
                    gBest = pValue;
                }

                ///存储画图数据
                positionArr[i][j][0] = particleArr[j].getX();
                positionArr[i][j][1] = particleArr[j].getY();
            }
        }

        //输出保留6位小数
        System.out.println("最优目标函数值:" + String.format("%.6f", gBest));
        System.out.println("最优x:" + String.format("%.6f", bestX));
        System.out.println("最优y:" + String.format("%.6f", bestY));
        System.out.println("求解时间:" + (System.currentTimeMillis() - start) + "ms");

        return positionArr;
    }

    /**
     * 目标函数
     *
     * @param x
     * @param y
     * @return
     */
    private double objectFunction(double x, double y) {
        //目标:在变量区间范围最小化 y=x^2+y^2-xy-10x-4y+60
        return Math.pow(x, 2) + Math.pow(y, 2) - x * y - 10 * x - 4 * y + 60;
    }

    /**
     * 粒子类
     */
    class Particle {
        private double x;
        private double y;
        //x,y坐标的上下限
        private double xMin;
        private double xMax;
        private double yMin;
        private double yMax;
        //x轴方向上的速度
        private double xV;
        //y轴方向上的速度
        private double yV;
        //该粒子找到的历史最优解
        private double pBest;
        //该粒子找到的历史最优解对应的x和y
        private double bestX;
        private double bestY;

        public Particle(double xMin, double xMax, double yMin, double yMax, double xV, double yV, Random random) {
            this.xMin = xMin;
            this.xMax = xMax;
            this.yMin = yMin;
            this.yMax = yMax;
            this.xV = xV;
            this.yV = yV;
            //初始化粒子信息
            this.initParticle(random);
        }

        /**
         * 初始化粒子信息
         * 即初始化位置
         */
        public void initParticle(Random random) {
            this.x = random.nextDouble() * (this.xMax - this.xMin) + this.xMin;
            this.y = random.nextDouble() * (this.yMax - this.yMin) + this.yMin;
//            System.out.println("this.x:" + this.x + "," + "this.y:" + this.y);
        }

        public double getX() {
            return x;
        }

        public void setX(double x) {
            this.x = x;
        }

        public double getY() {
            return y;
        }

        public void setY(double y) {
            this.y = y;
        }

        public double getxMin() {
            return xMin;
        }

        public void setxMin(double xMin) {
            this.xMin = xMin;
        }

        public double getxMax() {
            return xMax;
        }

        public void setxMax(double xMax) {
            this.xMax = xMax;
        }

        public double getyMin() {
            return yMin;
        }

        public void setyMin(double yMin) {
            this.yMin = yMin;
        }

        public double getyMax() {
            return yMax;
        }

        public void setyMax(double yMax) {
            this.yMax = yMax;
        }

        public double getxV() {
            return xV;
        }

        public void setxV(double xV) {
            this.xV = xV;
        }

        public double getyV() {
            return yV;
        }

        public void setyV(double yV) {
            this.yV = yV;
        }

        public double getpBest() {
            return pBest;
        }

        public void setpBest(double pBest) {
            this.pBest = pBest;
        }

        public double getBestX() {
            return bestX;
        }

        public void setBestX(double bestX) {
            this.bestX = bestX;
        }

        public double getBestY() {
            return bestY;
        }

        public void setBestY(double bestY) {
            this.bestY = bestY;
        }


    }

}

测试

package com.dam.heuristic.pso.test;

public class PsoMainRun {

    public static void main(String[] args) {
        PsoApi psoApi = new PsoApi(100, 2, 2, 3, 0.9, 1000);
        psoApi.solve();
    }

}

最优目标函数值:8.000000
最优x:8.000000
最优y:6.000000
求解时间:136ms

画图

package com.dam.heuristic.pso.test;

import com.dam.heuristic.vns.test.VnsApi;
import javafx.animation.KeyFrame;
import javafx.animation.Timeline;
import javafx.application.Application;
import javafx.geometry.Pos;
import javafx.scene.Scene;
import javafx.scene.canvas.Canvas;
import javafx.scene.canvas.GraphicsContext;
import javafx.scene.control.Button;
import javafx.scene.input.MouseEvent;
import javafx.scene.layout.BorderPane;
import javafx.scene.layout.HBox;
import javafx.scene.paint.Color;
import javafx.stage.Stage;
import javafx.util.Duration;

import java.io.File;
import java.io.FileInputStream;
import java.util.Arrays;

public class PsoPaint extends Application {

    //当前的时间轴
    private Timeline nowTimeline;
    //绘图位置坐标
    private double[][][] positionArr;

    public static void main(String[] args) {
        launch(args);
    }

    @Override
    public void start(Stage primaryStage) throws Exception {

        调用算法获取绘图数据
        PsoApi psoApi = new PsoApi(500, 2, 2, 3, 0.9, 6000);
        this.positionArr = psoApi.solve();

        画图
        try {
            BorderPane root = new BorderPane();
            root.setStyle("-fx-padding: 20;");
            Scene scene = new Scene(root, 1600, 900);
            double canvasWid = 800;
            double canvasHei = 800;
            //根据画布大小缩放坐标值
            this.fixPosition(canvasWid - 100, canvasHei - 100);

            //画布和画笔
            HBox canvasHbox = new HBox();
            Canvas canvas = new Canvas();
            canvas.setWidth(canvasWid);
            canvas.setHeight(canvasHei);
            canvasHbox.setPrefWidth(canvasWid);
            canvasHbox.getChildren().add(canvas);
            canvasHbox.setAlignment(Pos.CENTER);
            canvasHbox.setStyle("-fx-spacing: 20;" +
                    "-fx-background-color: #ecf1c3;");
            root.setTop(canvasHbox);
            GraphicsContext paintBrush = canvas.getGraphicsContext2D();

            //启动
            HBox hBox2 = new HBox();
            Button beginButton = new Button("启动粒子群仿真");
            hBox2.getChildren().add(beginButton);
            root.setBottom(hBox2);
            hBox2.setAlignment(Pos.CENTER);
            //启动仿真以及暂停仿真
            beginButton.addEventHandler(MouseEvent.MOUSE_CLICKED, event -> {
                nowTimeline.play();
            });

            //创建扫描线连接动画
            nowTimeline = new Timeline();
            createAnimation(paintBrush, 0.05);

            primaryStage.setScene(scene);
            primaryStage.show();
        } catch (Exception e) {
            e.printStackTrace();
        }
    }

    /**
     * 修正cityPositionArr的坐标,让画出来的点在画布内
     *
     * @param width
     * @param height
     */
    private void fixPosition(double width, double height) {
        double minX = Double.MAX_VALUE;
        double maxX = -Double.MAX_VALUE;
        double minY = Double.MAX_VALUE;
        double maxY = -Double.MAX_VALUE;

        for (int i = 0; i < this.positionArr.length; i++) {
            for (int j = 0; j < this.positionArr[0].length; j++) {
                minX = Math.min(minX, this.positionArr[i][j][0]);
                maxX = Math.max(maxX, this.positionArr[i][j][0]);
                minY = Math.min(minY, this.positionArr[i][j][1]);
                maxY = Math.max(maxY, this.positionArr[i][j][1]);
            }
        }

        double multiple = Math.max((maxX - minX) / width, (maxY - minY) / height);

        //转化为正数数
        for (int i = 0; i < this.positionArr.length; i++) {
            for (int j = 0; j < this.positionArr[0].length; j++) {
                if (minX < 0) {
                    this.positionArr[i][j][0] = this.positionArr[i][j][0] - minX;
                }
                if (minY < 0) {
                    this.positionArr[i][j][1] = this.positionArr[i][j][1] - minY;
                }
            }
        }

//        for (int i = 0; i < this.positionArr[0].length; i++) {
//            System.out.println(Arrays.toString(this.positionArr[99][i]));
//        }

        for (int i = 0; i < this.positionArr.length; i++) {
            for (int j = 0; j < this.positionArr[0].length; j++) {
                this.positionArr[i][j][0] = this.positionArr[i][j][0]/multiple;
                this.positionArr[i][j][1] = this.positionArr[i][j][1]/multiple;
            }
        }



    }

    /**
     * 用画笔在画布上画出所有的孔
     * 画第i代的所有粒子
     */
    private void drawAllCircle(GraphicsContext paintBrush, int i) {
        paintBrush.clearRect(0, 0, 2000, 2000);
        paintBrush.setFill(Color.RED);
        for (int j = 0; j < this.positionArr[i].length; j++) {
            drawCircle(paintBrush, i, j);
        }
    }

    /**
     * 用画笔在画布上画出一个孔
     * 画第i代的第j个粒子
     */
    private void drawCircle(GraphicsContext paintBrush, int i, int j) {
        double x = this.positionArr[i][j][0];
        double y = this.positionArr[i][j][1];
        double radius = 2;
        // 圆的直径
        double diameter = radius * 2;
        paintBrush.fillOval(x, y, diameter, diameter);
    }

    /**
     * 创建动画
     */
    private void createAnimation(GraphicsContext paintBrush, double speed) {
        for (int i = 0; i < this.positionArr[0].length; i++) {
            int finalI = i;
            KeyFrame keyFrame = new KeyFrame(Duration.seconds(i * speed), event -> drawAllCircle(paintBrush, finalI));
            nowTimeline.getKeyFrames().add(keyFrame);
        }
    }

}

粒子群的收敛过程可以参考下面的视频,刚开始时,粒子被随机放在x和y坐标都属于[-1000,1000]的任意位置,随着迭代次数的增加,粒子慢慢靠拢在一次,直到最后收敛于一个点。

在这里插入图片描述

Matlab作图

上面的图实在是太丑了,也是小编的JavaFx功底不咋行,下面改用Matlab进行绘图。

Java导出Excel数据

package com.dam.heuristic.pso.test;

import org.apache.poi.ss.usermodel.*;
import org.apache.poi.xssf.usermodel.XSSFCell;
import org.apache.poi.xssf.usermodel.XSSFCellStyle;
import org.apache.poi.xssf.usermodel.XSSFWorkbook;

import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.IOException;

public class PsoMainRun {

    public static void main(String[] args) {
        PsoApi psoApi = new PsoApi(500, 2, 2, 1, 0.5, 1000);
        double[][][] position = psoApi.solve();

        //创建WorkBook
        Workbook workbook = new XSSFWorkbook();
        Sheet xData = workbook.createSheet("xData");
        Sheet yData = workbook.createSheet("yData");

        //获取每一代,每一个粒子的x坐标
        CellStyle cellStyle = workbook.createCellStyle();
        //设置居中
        cellStyle.setAlignment(XSSFCellStyle.ALIGN_CENTER);
        for (int i = 0; i < position.length; i++) {
            Row xRow = xData.createRow(i);
            Row yRow = yData.createRow(i);
            for (int j = 0; j < position[0].length; j++) {
                Cell xCell = xRow.createCell(j);
                xCell.setCellValue(position[i][j][0]);
                //设置数据类型
                xCell.setCellType(XSSFCell.CELL_TYPE_NUMERIC);
                //设置居中
                xCell.setCellStyle(cellStyle);

                Cell yCell = yRow.createCell(j);
                yCell.setCellValue(position[i][j][1]);
                //设置数据类型
                yCell.setCellType(XSSFCell.CELL_TYPE_NUMERIC);
                //设置居中
                yCell.setCellStyle(cellStyle);

            }
        }

        try {
            FileOutputStream fileOutputStream = new FileOutputStream(new File("D:\\Desktop\\paintData.xlsx"));
            workbook.write(fileOutputStream);
            fileOutputStream.close();
            System.out.println("存储文件完成");
        } catch (FileNotFoundException e) {
            e.printStackTrace();
        } catch (IOException e) {
            e.printStackTrace();
        }

    }

}

在这里插入图片描述

数据样式
%% 读取数据
clear;
clc;
% 读取.mat数据
load xData
load yData

%% 绘制函数图像
figure 
% 确定画图区间
x = -1000:50:1000;
y = -1000:50:1000;
[x,y] = meshgrid(x,y);
z = x.^2 + y.^2 - x.*y - 10*x - 4*y + 60;
% 绘制网格
mesh(x,y,z)
% 加上坐标轴的标签
xlabel('x');  ylabel('y');  zlabel('z');  
% 冻结屏幕高宽比,使得一个三维对象的旋转不会改变坐标轴的刻度显示
axis vis3d 
 % 不关闭图形,继续在上面画粒子
hold on 

%% 绘制粒子图形
% 获取行数、列数
[r,c] = size(xData);
h=[];
for i = 1:r
    %获取每代所有粒子的x,y坐标
    xRow = xData(i,:);
    yRow = yData(i,:);
    zRow = xRow.^2 + yRow.^2 - xRow.*yRow - 10*xRow - 4*yRow + 60;
    if i==1
       % scatter3是绘制三维散点图的函数(这里返回h是为了得到图形的句柄,未来我们对其位置进行更新)
       h = scatter3(xRow,yRow,zRow,'*r');
       set(h,'XData',xRow,'YData',yRow,'ZData',zRow);
       pause(0.5);
    else
       %间隔0.04秒再画下一代
       pause(0.04);
       h.XData = xRow;
       h.YData = yRow;
       h.ZData = zRow; 
    end
end

在这里插入图片描述

  • 7
    点赞
  • 12
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
好的,我可以为您解答这个问题。 首先,让我们简要介绍一下粒子群算法(Particle Swarm Optimization,PSO)。 粒子群算法是一种进化计算技术,它模拟了鸟群或鱼群等生物群体的行为,通过不断地迭代来寻找最优解。每个个体在搜索空间中随机生成一组解,称为粒子。每个粒子都有一个位置和速度,位置表示当前解,速度表示搜索方向和速度大小。每个粒子还有一个适应度值,表示该解的优劣程度。 粒子群算法的主要思想是通过交换和合作来全局搜索解空间。每个粒子通过比较自己的适应度值和其周围粒子的适应度值来不断更新自己的位置和速度,直到找到最优解或达到最大迭代次数为止。 接下来,让我们来看一下如何使用粒子群算法求解 Rosenbrock 函数最小值。 Rosenbrock 函数是一个经典的非凸函数,其公式为: $f(x,y)=(1-x)^2+100(y-x^2)^2$ 该函数最小值为 $f(1,1)=0$。 我们可以将 Rosenbrock 函数的最小化问题转化为求解以下优化问题: $minimize\ f(x)=\sum_{i=1}^{n-1}[100(x_{i+1}-x_i^2)^2+(1-x_i)^2]$ 其中,$n$ 为 Rosenbrock 函数的维数。 下面是使用粒子群算法求解 Rosenbrock 函数最小值的 Python 代码示例: ```python import numpy as np # 定义 Rosenbrock 函数 def rosenbrock(x): return sum(100.0 * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0) # 粒子群算法 class PSO(): def __init__(self, dim, size, max_iter): self.w = 0.729 # 惯性权重 self.c1 = self.c2 = 1.49445 # 学习因子 self.dim = dim # 问题维度 self.size = size # 粒子数量 self.max_iter = max_iter # 最大迭代次数 self.X = np.random.uniform(-5, 5, (self.size, self.dim)) # 初始化粒子位置 self.V = np.random.uniform(-1, 1, (self.size, self.dim)) # 初始化粒子速度 self.pbest = self.X.copy() # 个体最优位置 self.gbest = self.X[np.argmin([rosenbrock(x) for x in self.X])] # 全局最优位置 def update(self): r1 = np.random.rand(self.size, self.dim) r2 = np.random.rand(self.size, self.dim) self.V = self.w * self.V + self.c1 * r1 * (self.pbest - self.X) + self.c2 * r2 * (self.gbest - self.X) self.X += self.V for i in range(self.size): if rosenbrock(self.X[i]) < rosenbrock(self.pbest[i]): self.pbest[i] = self.X[i] self.gbest = self.X[np.argmin([rosenbrock(x) for x in self.X])] def optimize(self): for i in range(self.max_iter): self.update() print('iter %d, best fitness: %f' % (i, rosenbrock(self.gbest))) pso = PSO(dim=2, size=50, max_iter=100) pso.optimize() print('global best:', pso.gbest) print('global best fitness:', rosenbrock(pso.gbest)) ``` 在上述代码中,我们首先定义了 Rosenbrock 函数,并实现了一个名为 PSO粒子群算法类。在类的初始化函数中,我们随机生成了初始粒子位置和速度,并将每个粒子的位置作为其个体最优位置。全局最优位置则初始化为粒子群中适应度值最小的粒子位置。 在算法的 update 函数中,我们按照粒子群算法的原理更新了粒子的速度和位置,并更新了每个粒子的个体最优位置和全局最优位置。在算法的 optimize 函数中,我们进行了多轮迭代,输出了每轮迭代的最优适应度值,最终输出了全局最优位置和最优适应度值。 当我们运行上述代码时,可以得到以下输出结果: ``` iter 0, best fitness: 7919.178681 iter 1, best fitness: 5281.189130 iter 2, best fitness: 3291.639364 ... iter 98, best fitness: 0.000000 iter 99, best fitness: 0.000000 global best: [1. 1.] global best fitness: 0.0 ``` 可以看到,经过多轮迭代,粒子群算法找到了 Rosenbrock 函数最小值 $f(1,1)=0$,并输出了全局最优位置和最优适应度值。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Hello Dam

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值