C#实现多目标粒子群优化算法（MOPSO）
这里如何用C#实现多目标粒子群优化算法。
本教程基于MATLAB实现多目标粒子群算法进行C#重现，MATLAB有许多现成的矩阵运算功能，C#需要自己写，因此该C#MOPSO源码在算法细节上不是很严谨，有很大优化改进空间，但程序运行没问题。
程序源码下载链接：
链接：https://pan.baidu.com/s/1PudUnIm_YHaMolNdGlzPxQ
提取码：hzdz
程序运行效果：
C#实现三维空间中点的绘制需要依赖一些外部库，如opengl等，因此本文直接将最终结果输出成txt格式，导入到MATLAB进行显示。
与MATLAB中一致的结果一致：
在有2个优化目标函数，并且优化目标函数设置合理的情况下，理想情况下，MOPSO的优化结果在平面内成线状。
在有3个优化目标函数，并且优化目标函数设置合理的情况下，理想情况下，MOPSO的优化结果在空间内成面状，如下图所示。
1、C#主程序
using System;
using System.IO;
struct empty_particle
{
public double[] Position;//记录粒子的位置
public double[] Velocity;//记录粒子的速度向量
public double[] Cost;//记录粒子的目标值向量
public double[] Best_Position;//粒子最佳位置向量
public double[] Best_Cost;//粒子最佳目标值向量
public int IsDominated;//粒子被支配个体向量
public int GridIndex;//粒子栅格索引向量
public int[] GridSubIndex;//粒子栅格子索引向量
};
struct empty_Grid
{
public double[] LB;//记录网络的经纬度
public double[] UB;//记录网络的经纬度
};
namespace MOPSO
{
class Program
{
static void Main(string[] args)
{
MOPSO();
}
//多目标粒子群优化
static void MOPSO()
{
int nVar = 2; //变量个数
int[] VarSize = { 1, nVar }; //变量矩阵大小
int VarMin = 0; //变量值定义域
int VarMax = 360; //注意:该函数变量不能出现负值
int MaxIt = 400; //最大迭代次数
int N = 40; //种群规模
int nRep = 50; //档案库大小
double w = 0.9; //惯性权重系数
double wdamp = 0.99; //惯性权重衰减率
double c1 = 1.7; //个体学习因子
double c2 = 1.8; //全局学习因子
int nGrid = 5; //每一维的分格数
double alpha = 0.1; //膨胀率
double beta = 2; //最佳选择压
double gamma = 2; //删除选择压
double mu = 0.1; //变异概率
empty_particle[] pop = new empty_particle[N]; // 初始粒子群
//初始化种群中的个体的位置、速度、目标值
for (int i = 0; i < N; i++)
{
//随机数
Random rd = new Random();
//初始化粒子的位置，即两个自变量
pop[i].Position = new double[] { rd.Next(VarMin*100, VarMax * 100)/100.0, rd.Next(VarMin * 100, VarMax * 100)/100.0 };
pop[i].Velocity = new double[] { 0.0, 0.0 };
pop[i].Cost = evaluate_objective(pop[i].Position);
pop[i].Best_Position = pop[i].Position;
pop[i].Best_Cost = pop[i].Cost;
}
//初始化种群中的个体的支配状况
DetermineDomination(pop);
//筛选出没有被支配的个体数
empty_particle[] rep = IsNotDominated(pop);
//创建初始网格
empty_Grid[] Grid = CreateGrid(rep, nGrid, alpha);
//找每个个体初始的网格位置
FindGridIndex(rep, Grid);
//%MOPSO主循环
for(int it = 0;it < MaxIt; it++)
{
//逐一个体更新速度和位置，0.5的概率发生变异
for (int i = 0; i<N; i++)
{
//从支配个体轮盘赌选出全局最佳个体
empty_particle leader = SelectLeader(rep, beta);
empty_particle[] Temp_rep = IsNotDominated(pop);
empty_particle[] new_rep = new empty_particle[rep.Length+ Temp_rep.Length];
rep.CopyTo(new_rep, 0);
Temp_rep.CopyTo(new_rep, rep.Length);
rep = new_rep;
Temp_rep = null;
new_rep = null;
Random rd = new Random();
double r1 = (rd.Next(0, 10000)) / 10000.0;
double r2 = (rd.Next(0, 10000)) / 10000.0;
// 速度更新
for (int j = 0;j< pop[i].Velocity.Length;j++)
{
pop[i].Velocity[j] = w * pop[i].Velocity[j] + c1 * r1 * (pop[i].Best_Position[j] - pop[i].Position[j]) + c2*r2* (leader.Position[j] - pop[i].Position[j]);
}
//位置更新
for (int j = 0; j < pop[i].Position.Length; j++)
{
pop[i].Position[j] = pop[i].Position[j] + pop[i].Velocity[j];
}
//限制变量变化范围
pop[i].Position = limitToPosition(pop[i].Position, VarMin, VarMax);
//计算目标函数值
pop[i].Cost = evaluate_objective(pop[i].Position);
//应用变异策略
double pm = (1 - (it - 0) / Math.Pow((MaxIt - 0),(1 / mu))); //变异概率逐渐变小
empty_particle NewSol;
NewSol.Position = Mutate(pop[i].Position, pm, VarMin, VarMax);
NewSol.Cost = evaluate_objective(NewSol.Position); // 计算变异后的目标值
if(Dominates(NewSol.Cost, pop[i].Cost))
{
pop[i].Position = NewSol.Position;
pop[i].Cost = NewSol.Cost;
}
else
{
double rand = (rd.Next(0, 10000)) / 10000.0;
if(rand < 0.5)
{
pop[i].Position = NewSol.Position;
pop[i].Cost = NewSol.Cost;
}
}
//% 如果当前个体优于先前最佳个体,则替换之
if(Dominates(pop[i].Cost, pop[i].Best_Cost))
{
pop[i].Best_Position = pop[i].Position;
pop[i].Best_Cost = pop[i].Cost;
}
else
{
double rand = (rd.Next(0, 10000)) / 10000.0;
if (rand < 0.5)
{
pop[i].Best_Position = pop[i].Position;
pop[i].Best_Cost = pop[i].Cost;
}
}
//Console.WriteLine(i);
}
DetermineDomination(rep);
//筛选出没有被支配的个体数
empty_particle[] T_rep = IsNotDominated(rep);
rep = T_rep;
Grid = CreateGrid(rep, nGrid, alpha);
//找每个个体初始的网格位置
FindGridIndex(rep, Grid);
if(rep.Length>nRep)
{
int Extra = rep.Length - nRep;
for(int e=0;e < Extra;e++)
{
rep = DeleteOneRepMemebr(rep, gamma);
}
}
Console.Write("迭代次数 =");
Console.WriteLine(it);
w = w * wdamp;
}
Console.WriteLine("迭代完成！！");
pc2TXT(rep);
}
//计算目标函数值
static double[] evaluate_objective(double[] x)
{
double f1 = x[0] * 0.01;
double f2 = (361 - x[0]) * (361 - x[1]) * 0.02;
double f3 = x[1] * 0.01;
double[] f = new double[] { f1, f2, f3 };
return f;
}
//判断全局支配状况,返回0 = 非支配解
static void DetermineDomination(empty_particle[] pop)
{
int nPop = pop.Length;
for(int i = 0;i<nPop;i++)
{
pop[i].IsDominated = 0;//初始化为互不支配
}
for(int i = 0;i< nPop - 1; i++)
{
for(int j = i+1;j< nPop;j++)
{
if(Dominates(pop[i].Cost, pop[j].Cost))
{
pop[j].IsDominated = 1;
}
if(Dominates(pop[j].Cost, pop[i].Cost))
{
pop[i].IsDominated = 1;
}
}
}
}
//%判断两个目标值x,y的支配状态
//x支配y,返回1;y支配x,返回0
static bool Dominates(double[] x,double[] y)
{
int flag_all = 1;
int flag_any = 0;
//x对应的每一位上的数要小于等于y
for (int i = 0;i<x.Length;i++)
{
if(x[i] > y[i])
{
flag_all = 0;
break;
}
}
//x对应的每一位上的数一定要比y中的某一位要小
for (int i = 0; i < x.Length; i++)
{
for (int j = 0; j < y.Length; j++)
{
if (x[i] < y[j])
{
flag_any++;
break;
}
}
}
if(flag_all == 1 && flag_any == x.Length)
{
return true;
}
else
{
return false;
}
}
//判断没有被支配的个体数
static empty_particle[] IsNotDominated(empty_particle[] pop)
{
int nPop = pop.Length;
int num = 0;
for (int i = 0; i < nPop; i++)
{
if (pop[i].IsDominated == 0)
num++;
}
empty_particle[] rep = new empty_particle[num];
int j=0;
for (int i = 0; i < nPop; i++)
{
if (pop[i].IsDominated == 0)
{
empty_particle temp = new empty_particle();
double temp_double;
int temp_int;
temp.Position = new double[pop[i].Position.Length];
for (int k=0;k<pop[i].Position.Length;k++)
{
temp_double = pop[i].Position[k];
temp.Position[k] = temp_double;
}
temp.Cost = new double[pop[i].Cost.Length];
for (int k = 0; k < pop[i].Cost.Length; k++)
{
temp_double = pop[i].Cost[k];
temp.Cost[k] = temp_double;
}
temp.Velocity = new double[pop[i].Velocity.Length];
for (int k = 0; k < pop[i].Velocity.Length; k++)
{
temp_double = pop[i].Velocity[k];
temp.Velocity[k] = temp_double;
}
temp.Best_Cost = new double[pop[i].Best_Cost.Length];
for (int k = 0; k < pop[i].Best_Cost.Length; k++)
{
temp_double = pop[i].Best_Cost[k];
temp.Best_Cost[k] = temp_double;
}
temp.Best_Position = new double[pop[i].Best_Position.Length];
for (int k = 0; k < pop[i].Best_Position.Length; k++)
{
temp_double = pop[i].Best_Position[k];
temp.Best_Position[k] = temp_double;
}
//temp.GridSubIndex = new int[pop[i].GridSubIndex.Length];
//for (int k = 0; k < pop[i].GridSubIndex.Length; k++)
//{
// temp_int = pop[i].GridSubIndex[k];
// temp.GridSubIndex[k] = temp_int;
//}
temp_int = pop[i].IsDominated;
temp.IsDominated = temp_int;
temp_int = pop[i].GridIndex;
temp.GridIndex = temp_int;
rep[j] = temp;
j++;
}
}
return rep;
}
//创建网格
static empty_Grid[] CreateGrid(empty_particle[] pop, int nGrid,double alpha)
{
int nPop = pop.Length;
int N = pop[0].Cost.Length;
//找出每个个体中Cost每一位上的最大最小值
double[] cmin = new double[3] { pop[0].Cost[0], pop[0].Cost[1], pop[0].Cost[2] };
double[] cmax = new double[3] { pop[0].Cost[0], pop[0].Cost[1], pop[0].Cost[2] };
for (int i = 0; i < nPop; i++)
{
for(int j = 0;j<3;j++)
{
if(pop[i].Cost[j]< cmin[j])
{
cmin[j] = pop[i].Cost[j];
}
if(pop[i].Cost[j] > cmax[j])
{
cmax[j] = pop[i].Cost[j];
}
}
}
double[] dc = new double[3] { cmax[0] - cmin[0], cmax[1] - cmin[1], cmax[2] - cmin[2] };
//扩大10%
for(int i=0;i<3;i++)
{
cmin[i] = cmin[i] - alpha * dc[i];
cmax[i] = cmax[i] + alpha * dc[i];
}
int nObj = cmin.Length;
//创建网格
empty_Grid[] Grid = new empty_Grid[3];
//大数
double big = 100000;
for(int i=0;i<nObj;i++)
{
double interval = (cmax[i] - cmin[i]) / nGrid;
Grid[i].LB = new double[nGrid+2];
Grid[i].LB[0] = -big;
for(int j =1;j< Grid[i].LB.Length;j++)
{
Grid[i].LB[j] = cmin[i] + interval * (j - 1);
}
Grid[i].UB = new double[nGrid + 2];
Grid[i].UB[Grid[i].LB.Length-1] = big;
for (int j = 0; j < Grid[i].LB.Length-1; j++)
{
Grid[i].UB[j] = cmin[i] + interval * j;
}
}
return Grid;
}
//判断所在的网格
static void FindGridIndex(empty_particle[] rep, empty_Grid[] Grid)
{
int nObj = rep[0].Cost.Length;
int nGrid = Grid[0].LB.Length;
for (int i = 0; i < rep.Length; i++)
{
rep[i].GridSubIndex = new int[nGrid];
for(int j=0;j<nObj;j++)
{
for(int k=0;k< nGrid;k++)
{
if(rep[i].Cost[j]<=Grid[j].UB[k])
{
rep[i].GridSubIndex[j] = k;
break;
}
}
}
rep[i].GridIndex = rep[i].GridSubIndex[0];
for(int j=1;j<nObj;j++)
{
rep[i].GridIndex = rep[i].GridIndex - 1;
rep[i].GridIndex = nGrid * rep[i].GridIndex;
rep[i].GridIndex = rep[i].GridIndex + rep[i].GridSubIndex[j];
}
}
}
//从全局支配个体中找出一个最佳个体
static empty_particle SelectLeader(empty_particle[] rep,double beta)
{
//获取所有个体的网格序号
int[] GI = new int[rep.Length];
for(int i=0;i< rep.Length;i++)
{
GI[i] = rep[i].GridIndex;
}
//对网格序号进行排序
int[] T_OC = GI;
Array.Sort(T_OC);
int T_length = 1;
//去除掉所有相同的序号
for (int i = 0; i < T_OC.Length - 1; i++)
{
if(T_OC[i]!=T_OC[i+1])
{
T_length++;
}
}
int[] OC = new int[T_length];
OC[0] = T_OC[0];
int T_i = 1;
for (int i = 0; i < T_OC.Length - 1; i++)
{
if (T_OC[i] != T_OC[i + 1])
{
OC[T_i] = T_OC[i+1];
T_i++;
}
}
//对相同的序号进行计数
int[] N = new int[T_length];
for (int i = 0; i < OC.Length; i++)
{
for(int j=0;j < T_OC.Length;j++)
{
if (OC[i] == T_OC[j])
{
N[i]++;
}
}
}
// 计算选择概率,为了增加多样性,尽量不选多次出现的个体
// 如果N大P就小, 即多次出现的栅格点被选中的概率小
double[] P = new double[T_length];
double T_sum = 0;
for (int i = 0; i < T_length; i++)
{
P[i] = Math.Pow(Math.E, -beta * N[i]);
T_sum = T_sum + P[i];
}
for (int i = 0; i < T_length; i++)
{
P[i] = P[i]/T_sum;
}
//轮盘赌策略选择
int sci = RouletteWheelSelection(P);
int sc = OC[sci]; //轮盘赌选择的栅格点
int T_SCM_Length = 0;
for (int i = 0; i < rep.Length; i++)
{
GI[i] = rep[i].GridIndex;
}
for (int i = 0; i < GI.Length; i++)
{
if (GI[i] == sc)
{
T_SCM_Length++;
}
}
int[] SCM = new int[T_SCM_Length];
int T_SCM_i = 0;
for(int i = 0; i < GI.Length; i++)
{
if (GI[i] == sc)
{
SCM[T_SCM_i] = i;
T_SCM_i++;
}
}
Random rd = new Random();
int smi = rd.Next(0, T_SCM_Length);
int sm = SCM[smi];
//当前全局最佳位置点
return rep[sm];
}
//%轮盘赌选择一个较好的支配个体
static int RouletteWheelSelection(double[] P)
{
//随机数
Random rd = new Random();
double r = (rd.Next(0,10000))/10000.0;
double sum = 0;
int i;
for(i=0;i<P.Length;i++)
{
sum = sum + P[i];
if(sum>r)
{
break;
}
}
if (i >= P.Length)
i = P.Length - 1;
return i;
}
//位置限制
static double[] limitToPosition(double[] Position, int VarMin,int VarMax)
{
for(int i =0;i< Position.Length;i++)
{
if (Position[i]< VarMin)
{
Position[i] = VarMin;
}
if(Position[i] > VarMax)
{
Position[i] = VarMax;
}
}
return Position;
}
//使用变异策略
static double[] Mutate(double[] x,double pm,int VarMin,int VarMax)
{
int nVar = x.Length;
Random rd = new Random();
int j = rd.Next(0,nVar - 1);
double dx = pm * (VarMax - VarMin);
double lb = x[j] - dx;
if(lb < VarMin)
{
lb = VarMin;
}
double ub = x[j] + dx;
if(ub > VarMax)
{
ub = VarMax;
}
double[] xnew = new double[x.Length];
for(int i=0;i< x.Length;i++)
{
double temp = x[i];
xnew[i] = temp;
}
xnew[j] = rd.Next((int)lb, (int)ub);
return xnew;
}
//删除档案库中的一个个体
static empty_particle[] DeleteOneRepMemebr(empty_particle[] rep, double gamma)
{
int[] GI = new int[rep.Length];
for(int i=0;i< rep.Length;i++)
{
GI[i] = rep[i].GridIndex;
}
//对网格序号进行排序
int[] T_OC = GI;
Array.Sort(T_OC);
int T_length = 1;
//去除掉所有相同的序号
for (int i = 0; i < T_OC.Length - 1; i++)
{
if (T_OC[i] != T_OC[i + 1])
{
T_length++;
}
}
int[] OC = new int[T_length];
OC[0] = T_OC[0];
int T_i = 1;
for (int i = 0; i < T_OC.Length - 1; i++)
{
if (T_OC[i] != T_OC[i + 1])
{
OC[T_i] = T_OC[i + 1];
T_i++;
}
}
//对相同的序号进行计数
int[] N = new int[T_length];
for (int i = 0; i < OC.Length; i++)
{
for (int j = 0; j < T_OC.Length; j++)
{
if (OC[i] == T_OC[j])
{
N[i]++;
}
}
}
// 计算选择概率,为了增加多样性,尽量不选多次出现的个体
// 如果N大P就小, 即多次出现的栅格点被选中的概率小
double[] P = new double[T_length];
double T_sum = 0;
for (int i = 0; i < T_length; i++)
{
P[i] = Math.Pow(Math.E,gamma * N[i]);
T_sum = T_sum + P[i];
}
for (int i = 0; i < T_length; i++)
{
P[i] = P[i] / T_sum;
}
//轮盘赌策略选择
int sci = RouletteWheelSelection(P);
int sc = OC[sci]; //轮盘赌选择的栅格点
int T_SCM_Length = 0;
for (int i = 0; i < rep.Length; i++)
{
GI[i] = rep[i].GridIndex;
}
for (int i = 0; i < GI.Length; i++)
{
if (GI[i] == sc)
{
T_SCM_Length++;
}
}
int[] SCM = new int[T_SCM_Length];
int T_SCM_i = 0;
for (int i = 0; i < GI.Length; i++)
{
if (GI[i] == sc)
{
SCM[T_SCM_i] = i;
T_SCM_i++;
}
}
Random rd = new Random();
int smi = rd.Next(0, T_SCM_Length);
int sm = SCM[smi];
empty_particle[] TT_rep = new empty_particle[rep.Length - 1];
for (int i = 0; i < sm; i++)
{
TT_rep[i] = rep[i];
}
for (int i = sm; i < TT_rep.Length; i++)
{
TT_rep[i] = rep[i + 1];
}
return TT_rep;
}
//将点云保存到txt
static void pc2TXT(empty_particle[] rep)
{
FileStream fs = new FileStream("D:\\12其他\\大师姐\\多目标优化\\C#\\MOPSO\\PC.txt", FileMode.Create);
StreamWriter sw = new StreamWriter(fs);
//开始写入
foreach (var item in rep)
{
sw.WriteLine(item.Cost[0] + "\t" + item.Cost[1] + "\t" + item.Cost[2]);
//清空缓冲区
sw.Flush();
}
//关闭流
sw.Close();
fs.Close();
}
}
}2、将结果输出到对应txt文件
将文件保存路径修改成自己的路径，并在主程序运行完后调用。
//将点云保存到txt
static void pc2TXT(empty_particle[] rep)
{
FileStream fs = new FileStream("D:\\12其他\\大师姐\\多目标优化\\C#\\MOPSO\\PC.txt", FileMode.Create);
StreamWriter sw = new StreamWriter(fs);
//开始写入
foreach (var item in rep)
{
sw.WriteLine(item.Cost[0] + "\t" + item.Cost[1] + "\t" + item.Cost[2]);
//清空缓冲区
sw.Flush();
}
//关闭流
sw.Close();
fs.Close();
}3、MATLAB中显示程序运行结果（MATLAB代码）
clc;
clear;
a = textread('PC.txt');%修改成自己的txt
figure(1);
location = a'; %取最优结果
scatter3(location(1,:),location(2,:),location(3,:),'filled','b');
xlabel('f1');ylabel('f2'); zlabel('f3');
title('Pareto 最优边界图');
box on;4、注意点！！！
如果C#跑出来的结果和MATLAB的结果不一致。
可能原因1：很有可能是C#中的随机函数并不随机导致的，因为C#中的随机取值是由时间决定的，因此有概率出现在循环中大量出现同一个随机数。
解决方案：(1)在循环的适当位置添加一些无用的输出，降低循环的速度。(对于部分机器没有用)；（2）把Random rd = new Random()拿出循环，设置rd为全局的变量。
可能原因2：gamma值设置过大，超出了double定义
解决方案：适当减小gamma值。
可能原因3：C#在数组传值时，不像MATLAB只传了数值，而是直接传了地址，导致在调整某个数组的时候，另一个数组也被修改。
解决方案：检查每个数组传值的地方是否有类似问题，并进行修改。
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
