基本人工蜂群算法的C#实现
用人工蜂群算法来实现对于高铁在区间上运行的能耗以及运行时间的计算
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
using System.Threading.Tasks;
using 小论文01_改进人工蜂群算法;
namespace _Solution
{
public class Solution
{
double[] fit_value = new double[100];//存放100个解向量的适应度
int maxIteration = 100;//最大的迭代个数为100代
double best_gap_time;//解向量中最佳的偏差时间
double run_time = 0;//平时线路上的运行时间
double makeup_time = 0;//坡度补偿时间
int best_solution_index = 0;//最佳偏差时间在矩阵中的行数
double[] best_velocity = new double[10001];//最佳的速度解向量
double[] failure = new double[100];//变异矩阵，若迭代后解的适应度未达到更好，则更新该矩阵，达到次数后即变异。
double Op_Energy = 0;//在速度解向量下列车运行的能量消耗
double Po_Energy = 0;//列车在区段运行时，重力势能的变化量
public double[] extractor(double[,] v, int row)//将矩阵的一行（一个解向量）提取出来，方便后续迭代操作。
{
double[] result = new double[10001];
for (int i = 0; i < 10001; i++)
{
result[i] = v[row, i];
}
return result;
}
public int find_best(double[] v)//在数组中找最小值，并返回最小值所在下标
{
double min = v[0];
int index = 0;
for (int i = 0; i < v.Length; i++)
{
if (v[i] < min)
{
min = v[i];
index = i;
}
}
return index;
}
public int randomIndex_generator()//在0-100的范围内随意生成下标
{
Random rand = new Random();
int index = rand.Next(0, 100);
return index;
}
static void Main()
{
Solution s = new Solution();
Calculation cal = new Calculation();
double[] temp = new double[10001];
double[,] solution_matrix = new double[100, 10001];
double[] time = new double[100];
Random r = new Random();
Random rd = new Random();
Fit f = new Fit();
double[] fit_value = new double[100];
//先用0将二维数组填满
for (int i = 0; i < 100; i++)
{
for (int j = 0; j < 10001; j++)
{
solution_matrix[i, j] = 0;
}
}
//再用最大加速度填满加速阶段和减速阶段
for (int i = 0; i < 100; i++)
{
for (int j = 1; j < 26; j++)
{
solution_matrix[i, j] = Math.Sqrt(solution_matrix[i, j - 1] * solution_matrix[i, j - 1] + 2 * cal.a_lim * 12.493);
}
for (int k = 1; k < 26; k++)
{
solution_matrix[i, k] *= 3.6;
}
}
for (int i = 0; i < 100; i++)
{
for (int j = 9999; j > 9975; j--)
{
solution_matrix[i, j] = Math.Sqrt(2 * cal.a_lim * 12.493 + solution_matrix[i, j + 1] * solution_matrix[i, j + 1]);
}
for (int k = 9999; k > 9975; k--)
{
solution_matrix[i, k] *= 3.6;
}
}
//最后将矩阵中的剩余位置用220-300的范围内的随机数添满
for (int i = 0; i < 100; i++)
{
for (int j = 26; j < 9975; j++)
{
solution_matrix[i, j] = 220 + r.NextDouble() * (cal.V_max - 220);
}
}
//根据生成的100组初始解，计算出每组的偏差时间，并放入偏差时间的数组中储存
for (int i = 0; i < 100; i++)
{
for (int j = 0; j < 10001; j++)
{
temp[j] = solution_matrix[i, j];
}
time[i] = cal.gap_time(temp);
}
//计算每一组初始解的适应度函数值，并放入适应度数组中保存
for (int i = 0; i < 100; i++)
{
for (int j = 0; j < 10001; j++)
{
temp[j] = solution_matrix[i, j];
}
fit_value[i] = f.fitValue(temp);
}
//将舍弃旧解的标记数组全部初始为0
for (int i = 0; i < 100; i++)
{
s.failure[i] = 0;
}
//开始迭代
int row, column;
for (int i = 0; i < s.maxIteration; i++)//开始迭代
{
for (row = 0; row < 100; row++)
{
temp = s.extractor(solution_matrix, row);//将解矩阵中的一个解向量取出，放到临时的数组中保存
for (column = 26; column < 9976; column++)
{
int pos = s.randomIndex_generator();//随机生成一个下标；
double tNumber = temp[column] + r.Next(-1, 2) * (Math.Abs(temp[column] - solution_matrix[row, pos]));//在初始解的邻域附近找新的可行解
temp[column] = tNumber;
if (fit_value[row] < f.fitValue(temp))//如果新解的适应度比原有解大，则接受原有的解，并更新变异矩阵
{
temp[column] = solution_matrix[row, column];
s.failure[row] += 1;
}
else
{
s.failure[row] = 0;
}//否则的话则让变异矩阵回归到0；且接受新的可行解
if (s.failure[row] > 6000)//当变异矩阵数值大于6000时（替换6000个分量后，仍没有提升）则全新生成一个新的解，替代原来的解，并更新适应度矩阵
{
for (int j = 26; j < 9976; j++)
{
temp[j] = 220 + r.NextDouble() * (cal.V_max - 220);
}
fit_value[row] = f.fitValue(temp);
for (int k = 0; k < 10001; k++)
{
solution_matrix[row, k] = temp[k];
}
break;
}
for (int m = 0; m < 10001; m++)//将新的解向量放回到矩阵中，并更新
{
solution_matrix[row, m] = temp[m];
}
fit_value[row] = f.fitValue(temp);//更新适应度值矩阵
time[row] = cal.gap_time(temp);//更新偏差时间矩阵
}
s.best_solution_index = s.find_best(fit_value);//找出最优的解所在的下标
s.best_gap_time = time[s.best_solution_index];//更新全局最优的时间
for (int q = 0; q < 10001; q++)//一次迭代完成后，更新最优的速度解向量
{
s.best_velocity[q] = solution_matrix[s.best_solution_index, q];
}
s.run_time = cal.run_time(s.best_velocity);
s.Op_Energy = cal.Operation_Energy(s.best_velocity);//计算最优速度解向量下，列车运行所消耗的能量
s.Po_Energy = cal.potentialEnergy();//计算区间重力势能变化量
s.makeup_time = (s.Po_Energy / (s.Po_Energy + s.Op_Energy)) * cal.plan_time;
//if (s.fit_value[s.best_solution_index] == Double.MinValue)
//{
// break;
//}
using (StreamWriter sw = new StreamWriter(@"D:\result\" + $"{i}new.txt"))//将最优的速度解向量输出到txt文本
{
sw.WriteLine("平时线路上运行消耗的能量为：" + s.Op_Energy);
sw.WriteLine("区间的重力势能变化量为：" + s.Po_Energy);
sw.WriteLine("平直线路上的运行时间为：" + s.run_time);
sw.WriteLine("坡度补偿时间为：" + s.makeup_time);
sw.WriteLine();
for (int p = 0; p < 10001; p++)
{
sw.WriteLine(s.best_velocity[p]);
}
}
}
}
Console.WriteLine("计算完成！");
}
}
}
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace 小论文01_改进人工蜂群算法
{
public class Calculation
{
public double M = 380000;//列车质量：满载时380000kg（单位：kg）
public double g = 9.8;//重力加速度9.8
public double S = 124930;//线路全长124930m（单位：m）
public double n = 10000;//分成10000个小区间
public double V_max = 300;//区间上的最大限速为300km/h
public double V_min = 0;//最小速度为0
public double a_lim = 0.6 * 9.8;//最大加（减）速度取人类所感觉舒适的临界加速度，这里取0.6g
public double plan_time = 34 * 60;//图定运行时间34min（单位：s）
public double[] slope = new double[] { 0.0, 3.5, 5.5, -2, 7, -10, -3.5, 6, -5, 4.5, 3.5, -6, 5, -4, -1, 0.0, 4.5, 6, 3.5, -5.5, 3, -4, 3.5, 3, -2, -1.5, 2.5, 3.5, 0.0, -4.5, -3, 3, -5.5, -6, 1.5, -3, -4.5, 5, 2, -1.5, -3.5, 4.5, -5, 10, -5.5, -3.5, 4.5, 3.5, -4, 0.0 };//线路上的具体坡度
public double[] distance = new double[] { 2350, 2420, 1780, 3120, 2570, 2420, 2170, 3230, 1760, 2450, 3560, 2010, 2300, 1890, 1710, 2300, 2640, 3150, 2490, 1970, 2130, 2980, 1760, 2430, 3370, 1250, 2360, 2780, 3060, 2090, 2430, 2760, 2730, 2730, 3140, 2960, 1780, 2450, 3430, 3250, 2960, 2870, 2730, 3030, 1950, 2840, 3120, 2760, 2840, 2530 };//线路的具体纵断面长度
public double kineticEnergy(double v)//动能的计算公式
{
double result = 0.5 * 380000 * (v/3.6) * (v/3.6);
return result;
}
public double primitiveFunc(double v)//总能量中除去动能外，其他部分的计算公式
{
double result = 0.5 * 2457.84 * (v/3.6) * (v/3.6) + (1 / 3) * 9.1238 * (v/3.6) * (v/3.6) * (v/3.6) + 0.25 * 0.491568 * v * v * v * v;
return result;
}
public double accelerateEnergy(double vs, double ve)//加速区间的能量消耗,vs为初速度；ve为末速度
{
double result, accleration;
accleration = ((ve/3.6) * (ve/3.6) - (vs/3.6) * (vs/3.6)) / (2 * 12.493);
result = (kineticEnergy(ve) - kineticEnergy(vs)) + (1 / accleration) * (primitiveFunc(ve) - primitiveFunc(vs));
return result;
}
public double decelerateEnergy(double vs, double ve)//减速区间的能量消耗，vs为初速度；ve为末速度，由于做功为负，但是所求能量消耗为正
{
double result, deceleration;
deceleration = ((vs/3.6) * (vs/3.6) - (ve/3.6 )* (ve/3.6)) / (2 * 12.493);
result = Math.Abs((kineticEnergy(ve) - kineticEnergy(vs)) - (1 / deceleration) * (primitiveFunc(ve) - primitiveFunc(vs)));
return result;
}
public double uniformEnergy(double v)//计算匀速区间所消耗能量
{
double result;
result = (2457.84 + 9.1238 * v + 0.491568 * v * v) * (124930 / 10000);
return result;
}
public double potentialEnergy()//由线路纵断面计算出的总的重力势能
{
double result, sum = 0;
for (int i = 0; i < slope.Length; i++)
{
double temp = slope[i] / 1000;
double r = temp / (Math.Sqrt(temp * temp + 1));
sum = sum + r * distance[i];
}
result = M * g * sum;
return result;
}
public double run_time(double[] v)//计算平直线路上的运行速度
{
double temp = 0, sum= 0, result = 0;
for (int i = 0; i < v.Length - 1; i++)
{
temp = (2 * 12.493) / ((v[i] / 3.6) + (v[i + 1] / 3.6));
sum = sum + temp;
}
result = sum / 60;
return result;
}
public double gap_time(double[] v)//计算偏差时间的函数，返回值单位：分钟
{
double temp = 0, run_time = 0, sum_operationEnergy = 0, Ep = 0, makeup_time, result = 0;
for (int i = 0; i < v.Length - 1; i++)
{
temp = (2 * (S / 10000)) / ((v[i] / 3.6) + (v[i + 1] / 3.6));
run_time = run_time + temp;
if (v[i] < v[i + 1])
{
sum_operationEnergy = sum_operationEnergy + accelerateEnergy(v[i] / 3.6, v[i + 1] / 3.6);
}
else if (v[i] > v[i + 1])
{
sum_operationEnergy = sum_operationEnergy + decelerateEnergy(v[i] / 3.6, v[i + 1] / 3.6);
}
else
{
sum_operationEnergy = sum_operationEnergy + uniformEnergy(v[i] / 3.6);
}
}
Ep = potentialEnergy();
makeup_time = (Ep / (Ep + sum_operationEnergy)) * plan_time;
result = Math.Abs(plan_time - (makeup_time + run_time)) / 60;
return result;
}
public double Operation_Energy(double[] v)//计算平直线路上列车运行所消耗的能量
{
double sum = 0;
for (int i = 0; i < v.Length - 1; i++)
{
if (v[i] < v[i + 1])
{
sum = sum + accelerateEnergy(v[i] / 3.6, v[i + 1] / 3.6);
}
else if (v[i] > v[i + 1])
{
sum = sum + decelerateEnergy(v[i] / 3.6, v[i + 1] / 3.6);
}
else
{
sum = sum + uniformEnergy(v[i] / 3.6);
}
}
return sum;
}
public double[] section_operationTime(double[] v)//计算每个小区间的运行时间
{
double[] result = new double[10000];
for (int i = 0; i < v.Length - 1; i++)
{
result[i] = 2 * (S / 10000) / (v[i] + v[i + 1]);
}
return result;
}
}
}
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace 小论文01_改进人工蜂群算法
{
public class Fit
{
public double punishValue(double[] v)//计算惩罚值
{
int count;
double L1, L2,L3;
double temp;
double result=0;
Calculation Cal = new Calculation();
for (int i = 0; i < v.Length-1; i++)
{
count = 0;L1 = L2 = L3 = 0;//让参数归0
if (v[i] > Cal.V_max)//判断是否违反了最大速度的约束
{
count += 1;
L1 = v[i] - Cal.V_max;
}
if (v[i] < Cal.V_min)//判断是否违反了最小速度的约束
{
count += 1;
L2 = Cal.V_min - v[i];
}
temp = (2 * (Cal.S / 10000)) / (v[i] + v[i + 1]);//计算每个小区间上的运行时间
if ((Math.Abs(v[i] - v[i + 1]) / temp) > Cal.a_lim)//判断在每一个区间上是否违反了加速度约束
{
count += 1;
L3= (Math.Abs(v[i] - v[i + 1]) / temp)-Cal.a_lim;
}
result = result+count * 100 + 100 * L1 + 100 * L2 + 100 * L3;//对解向量的每一个分量都进行判断，每违反一个约束都进行惩罚，且违反的越多惩罚的越多
}
return result;
}
public double fitValue(double[] v)
{
double result, temp, f;
Calculation Cal = new Calculation();
f = Cal.gap_time(v);
temp = punishValue(v);
result = temp + f;
return result;
}
}
}
扫一扫
2085
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
