纯C，使用遗传算法对病态线性方程组求解AX=B

原创

： 纯C，使用遗传算法对病态线性方程组求解AX=B

纯C，使用遗传算法对病态线性方程组求解AX=B

使用遗传算法对病态线性方程组求解

对于什么是遗传算法，在松鼠科学会上有一篇有趣的文章
遗传算法构建firefox图标

遗传算法 的大致框架如下图所示
下面的长代码有什么用呢？
它使用算法遗传解决了

AX=B

的问题

我就觉得效率不是很高，使用起来有很大的不便。

什么是遗传算法呢？

简单地说，遗传算法是一种解决问题的方法。它模拟大自然中种群在选择压力下的演化，从而得到问题的一个近似解。

将以下数据放入A.txt,B.txt
将上述文件放入程序文件同目录

A的数据

(1, 1, 1)
(2, 8, -1)
(2, 3, -1)
(2, 2, 4)
(3, 9, -1)
(3, 2, -1)
(3, 3, 4)
(4, 4, 1)
(5, 2, -1)
(5, 6, -1)
(5, 5, 4)
(6, 2, -1)
(6, 5, -1)
(6, 6, 4)
(7, 7, 1)
(8, 6, -1)
(8, 9, -1)
(8, 8, 4)
(9, 5, -1)
(9, 8, -1)
(9, 9, 4)
(10, 10, 1)
(11, 7, -1)
(11, 12, -1)
(11, 11, 4)
(12, 2, -1)
(12, 11, -1)
(12, 12, 4)

B的数据

(1, 1, 2)
(2, 1, 5)
(3, 1, 4)
(4, 1, 2)
(5, 1, 5)
(6, 1, 4)
(7, 1, 2)
(8, 1, 7)
(9, 1, 3)
(10, 1, 4)
(11, 1, 0)
(12, 1, 4）

X的结果

2.039260
2.469905
2.118222
2.037000
2.412758
2.226262
2.044251
2.806577
2.044782
4.032598
1.827159
5.236275

使用以下代码前需要把第6行和第10行的opencv库调用删除，如果没有把opencv的库装到某些奇怪的地方的话。
此代码共计600行阅读可能存在一定的困难

#define _CRT_SECURE_NO_WARNINGS
#include<stdio.h>
#include<math.h>
#include<string.h>
#include<stdlib.h>
#include<opencv.hpp>
#include <iostream>
#include <time.h>
#define _CRT_SECURE_NO_WARNINGS
using namespace cv;
#define SUM 200   //总共的染色体数量
#define MAXloop  2000000//最大循环次数
#define error 10   //若两次最优值之差小于此数则认为结果没有改变
#define crossp 0.1      //交叉概率
#define mp 0.4           //变异概率
#define length 12
float X[length+1];
//#define d 0
float a[length+1];
typedef struct NODE//稀疏方程组数据
{
	int row;
	int col;
	int date;
	struct NODE* next;
}node;
//   Ax = b    x的各元素在（0,255）之间的最小值 
struct gen                //定义染色体结构X
{
	float  info[length+1];        		//染色体结构 //表示length 个(0,255)的整数
	float  suitability;		//次染色体所对应的适应度函数值，在本题中为f(x)=|AX-B|的无穷范数 理想状态应该是1 选择时最大的  我们把它设置成整数就是f(x)的值
};
struct gen gen_group[SUM];//定义一个含有20个染色体的组
struct gen gen_new[SUM];
struct gen gen_yi[SUM];
struct gen gen_result;    //记录最优的染色体
int result_unchange_time=0; //记录在error前提下最优值为改变的循环次数
struct log                //形成链表，记录每次循环所产生的最优的适应度
{
	int suitability;
	struct log* next;
}llog, * head, * end;
int randNext(int left, int right);
int randsign(float p)//按概率p返回1
{
	if (randNext(0,32768) > (p * 32768))//32768是INT的最大值 
		return 0;
	else return 1;
}
int log_num;              //链表长度
node * head1, * head2;
void initiate();          	//初始化函数，主要负责产生初始化种群 
void evaluation(int flag);	//评估种群中各染色体的适应度，并据此进行排序 
void cross();				//交叉函数 
void selection(int x);			//选择函数 
int  record();				//记录每次循环产生的最优解并判断是否终止循环 
void mutation();			//变异函数 
float f(gen G);

void initiate()//初始化 
{
	int i;
	int j;
	for (i = 0; i < SUM; i++)
	{
		for (j = 1; j <= length; j++) {
				gen_group[i].info[j] =1.0* randNext(0,500000)/10000.0;//从1开始
		}
	}
	gen_result.suitability = 10;//最优解.对应的函数值 
}
//返回[left, right]的随机数
int randNext(int left, int right)
{
	static unsigned int seed = 0;
	seed++;
	srand((unsigned)time(NULL) + seed * seed);
	return rand() % (right - left + 1) + left;
}
int ff(int x,int y) {//防止近亲交配
	int count = 0;
	int i, j;
	for (i = 1;i <= length;i++) {
		if (1.0*gen_group[x].info[i] - 1.0*gen_group[y].info[i]<=0.1&& 1.0*gen_group[x].info[i] - 1.0*gen_group[y].info[i] >=- 0.1) {
			count++;
		}
	}
	if (count == length) {
		return 1;
	}
	else {
		return 0;
	}
}
void cross()	//交叉函数 交叉变异 交叉公式为x3=A * x1 +(1-A) * x2
{
	int i, j, k;
	int c;
	float A;
	//int mask1, mask2;
	int a[SUM];//sum=20 
	//printf("55");
	for (i = 0; i < SUM; i++)  a[i] = 0;//记录是否交叉过 0-19
	k = 0;		
	int x = 0;
	int s;
	for (i = 0; i < SUM; i++)
	{
		if (a[i] == 0)
		{
			j = randNext(i + 1, SUM - 1);
				while (ff(i,j)||a[j]!=0) {
					j = randNext(i + 1, SUM - 1);//产生一个在i+1，sum-1两个数之间的随机整数
				}
				//if (a[j] == 0)	break; //直到选择到未交叉过的
				A = (1.0*gen_group[i].suitability/(gen_group[i].suitability*1.0+1.0*gen_group[j].suitability));
				x = int(A * 1.0*length);
			if (randsign(crossp) == 1)		//按照crossp的概率对选择的染色体进行交叉操作 返回0或1   其值为1的概率为crossp  算数交换
			{

				for (c =1;c <= length;c++) {//交叉中
					gen_new[k].info[c] = ((1.0-A)*gen_group[i].info[c]+(1.0*A)*gen_group[j].info[c]);
					gen_new[k+1].info[c] = ((1.0-A) * gen_group[j].info[c] + ( 1.0*A) * gen_group[i].info[c]);
				}
				
				k = k + 2;
			}
			else 		//物理交换
			{
									for (c = 1;c <=x;c++) {//交叉中
										gen_new[k].info[c] = gen_group[i].info[c];
										gen_new[k + 1].info[c] = gen_group[j].info[c];
									}
									for (c = x + 1;c <= length;c++) {
										gen_new[k].info[c] = gen_group[j].info[c];
										gen_new[k + 1].info[c] = gen_group[i].info[c];
									}
				k = k + 2;
			}
			a[i] = a[j] = 1;//进行过交叉 
		}
	}
}
float  f(gen G) {         //计算||AX-B||的无穷范数   提供gen[x].    B的数值   A的所有数值
	 int i, j;
	float  sum = 0;
	float  max = 0;
	i = j = 1;
	node* p, * q;
	p = head1;
	//printf("19");
	q = head2;
	while (q != NULL&&i<length) {
			while (p->row == i&&p!=NULL) {		//按行进行计算  第i列对应x的第p->col行    乘积相加  再将结果减去B对应的值   q->date  //i<=196608
					sum =sum+1.0*p->date * 1.0*G.info[p->col]*1.0;
					p = p->next;
			}
			i++;//下一行
			sum -= (q->date)*1.0;//AX-B b的值是依次排列的
			if (sum < 0) {//绝对值
				sum = -sum;
			}
			//printf("19");
			q = q->next;
			if (sum > max) {
				max = sum;//max的理想状态应该是0
			}
			sum=0; //重置sum计算下一个
			if (p == NULL||q==NULL) {
				break;
			}
	}
	return max;
}
void evaluation(int flag)//评估种群中各染色体的适应度，并据此进行排序 
{
	int i, j;
	printf("333\n");
	struct gen* genp;
	struct gen * genp2;
	genp2 = (struct gen*)malloc(sizeof(gen));
	//float gentsuitability;
	//float x;
	int x;
	float k;
	if (flag == 0)			// flag=0的时候对父种群进行操作 ///genp = gen_group;
	{
		for (i = 0; i < SUM; i++)//计算各染色体对应的表达式值
		{
			gen_group[i].suitability = 10000.0 / (1.0 + f(gen_group[i]));     //在f最小  f=0时取最大值 大到小排序  f越小越接近理想值 即适应度越高
		}
		printf("999");
		for (i = 0; i < SUM - 1; i++)//按表达式的值进行排序，
		{
			for (j = i + 1; j < SUM; j++)
			{
				if (gen_group[i].suitability < gen_group[j].suitability)//对应的函数值进行比较 交换从大到小冒泡排序
				{
					for (int i1 = 1;i1 <= length;i1++) {
						k = gen_group[i].info[i1];
						gen_group[i].info[i1] = gen_group[j].info[i1];
						gen_group[j].info[i1] = k;
					}
					k = gen_group[i].suitability;
					gen_group[i].suitability = gen_group[j].suitability;
					gen_group[j].suitability = k;
				}
			}
		}
	}
	else if(flag==1)//genp = gen_new;
	{
			for (i = 0; i < SUM; i++)//计算各染色体对应的表达式值
			{
				gen_new[i].suitability =10000.0/(1.0+ f(gen_new[i]));     //在f最小  f=0时取最大值 大到小排序
			} 
		//	printf("9");
			for (i = 0; i < SUM - 1; i++)//按表达式的值进行排序，
			{
				for (j = i + 1; j < SUM; j++)
				{
					if (gen_new[i].suitability < gen_new[j].suitability)//对应的函数值进行比较 交换从大到小冒泡排序 
					{
						for (int i1 = 1;i1 <= length;i1++) {
							k = gen_new[i].info[i1];
							gen_new[i].info[i1] = gen_new[j].info[i1];
							gen_new[j].info[i1] = k;
						}
						k = gen_new[i].suitability;
						gen_new[i].suitability = gen_new[j].suitability;
						gen_new[j].suitability = k;
					}
				}
			}
	}
	free(genp2);
}
void selection(int x)//选择函数 选择以后的种群是有序
{
	int i, j, k;
	j = 0;
	i = SUM / 2 - 1;	//i=9 
	printf("11");
	if (x == 1)
	{
		for (i = 0;i <= SUM/2-5;i++) {
			for (j = 1;j <= length;j++)
			{
				gen_group[SUM/2-1+ i].info[j] = gen_yi[i*2].info[j];
			} 
			gen_group[SUM/2-1+ i].suitability = gen_yi[i*2].suitability;
		}
	}
	else if (x == 0) //genp = gen_new;
		{
			
		for (i = 0;i <=SUM/4;i++) {
			for (j = 1;j <= length;j++)
			{
				gen_group[i+SUM/4].info[j] = gen_new[i].info[j];
			}
			gen_group[i+SUM/4].suitability = gen_new[i].suitability;
		}
		}
}
int record()	//记录最优解和判断是否满足条件      14.662757
{
	float x;
	if (gen_group[0].suitability- 10000.0<=100&& gen_group[0].suitability - 10000.0 >=-100)
		return 1;
	x = gen_result.suitability - gen_group[0].suitability;//最优的染色体和gen0的差值 
	if (x < 0)	x = -x;
	if (x < error)//若两次最优值之差小于此数则认为结果没有改变
	{
		result_unchange_time++;//没有改变 
		if (result_unchange_time >= 100000)
			return 1;//连续100次无改变 满足条件 
	}
	else
	{
		gen_result = gen_group[0];
		gen_result.suitability = gen_group[0].suitability;
		result_unchange_time = 0;
	}
	return 0;
}
void mutation()//变异函数 变异概率mp=0.04
{
	int i, j;
	int k = 1;
	float gmp;
	float k1=0;
	struct gen* gen2 = (gen*)malloc(sizeof(gen));
	//gmp = 1 - pow(1 - mp, 11);//在基因变异概率为mp时整条染色体的变异概率  mp 0.04           //变异概率 gmp=1-(0.96)^11
	for (i = 0;i <= SUM - 1;i++) {
		for (j = 1;j <= length;j++) {
			gen_yi[i].info[j] = gen_group[i].info[j];
		}
	}
	for (i = 0; i < SUM; i++)
	{
		for (j = 0;j <=length/3;j++) 
		{
					k = randNext(1, length);
					gen_yi[i].info[k] = 1.0*randNext(0, 50000)/1000.0;//变异 
		}
		 gen_yi[i].suitability = 10000.0/(1.0+1.0*f(gen_yi[i]));
	}
	for (i = 0; i < SUM - 1; i++)//重新排序
	{
		for (j = i + 1; j < SUM; j++)
		{
			if (gen_yi[i].suitability < gen_yi[j].suitability)
			{
				for (int i1 = 1;i1 <= length;i1++) {
					k1 = gen_yi[i].info[i1];
					gen_yi[i].info[i1] = gen_yi[j].info[i1];
					gen_yi[j].info[i1] = k1;
				}
				 k1= gen_yi[i].suitability;
				gen_yi[i].suitability = gen_yi[j].suitability;
				gen_yi[j].suitability = k1；
			}
		}
	}
	//printf("12\n");
}
node* createList(char str[]);//读数据 
void destroyList(node* head);
void suanfa(node* head1, node* head2);
//Mat image(Size(4, 1), CV_8UC3,Scalar(255, 255, 255));
void sort() {
	int i, j;
	int x;
	float k;
	gen* genp2 = (struct gen*)malloc(sizeof(gen));
	for (i = 0; i < SUM - 1; i++)//按表达式的值进行排序，
	{
		for (j = i + 1; j < SUM; j++)
		{
			if (gen_group[i].suitability <gen_group[j].suitability)//对应的函数值进行比较 交换从大到小冒泡排序 
			{
				for (int i1 = 1;i1 <= length;i1++) {
					k = gen_group[i].info[i1];
					gen_group[i].info[i1] = gen_group[j].info[i1];
					gen_group[j].info[i1] = k;

				}
				k = gen_group[i].suitability;
				gen_group[i].suitability = gen_group[j].suitability;
				gen_group[j].suitability = k;
			}
		}
	}
	free(genp2);
}
int main()
{
	int i, flag;
	char str1[] = "A.txt";
	char str2[] = "B.txt";
	//printf("1");
	head1 = createList(str1);//A
	head2 = createList(str2);//B
	//printf("19");
	//以下是遗传算法的内容
	{
				initiate();			//产生初始化种群 
				printf("19");
				evaluation(0);		//对初始化种群进行评估、排序
				printf("15");
				for (i = 0; i < MAXloop; i++)//MAXloop 1200
			    {
					cross();			//进行交叉操作 产生new
					printf("19");
					evaluation(1);	//对子种群进行评估、排序 对new排序
					selection(0);		//对父子种群中选择最优的NUM个作为新的父种群 
					sort();
					printf("19");
					if (record() == 1)	//满足终止规则1，则flag=1并停止循环 
					{
						flag = 1;
						break;
					}
					printf("19");
					mutation();			//变异操作  自带评估排序功能
					selection(1);		//对父子种群中选择最优的NUM个作为新的父种群
					sort();
					//if (i % 10 == 0)
					{
						printf("\n number = %d  \n", i);
						printf(" \ngengroup  %f  %f  %f\n", gen_group[0].suitability,gen_group[SUM/2].suitability,gen_group[SUM-1].suitability);
						printf(" \ngenyi  %f  %f  %f\n", gen_yi[0].suitability,gen_yi[SUM/2].suitability, gen_yi[SUM - 1].suitability);
						printf(" \ngennew  %f   %f %f\n", gen_new[0].suitability,gen_new[SUM/2].suitability, gen_new[SUM - 1].suitability);
						printf("\n gen_result %f  \n", gen_result.suitability);
					}
				}
	}
	suanfa(head1, head2);
	destroyList(head1);
	destroyList(head2);
	//imwrite("X.jpg", image);
	return 0;
}
node* createList(char str[])//读文件
{
	FILE* fp;
	fp = fopen(str, "r");
	node* p, * q, * head;
	head = (node*)malloc(sizeof(node));
	q = head;
	int rows, cols, udate;
	printf("%c", fgetc(fp));
	while (!feof(fp))
	{
		fscanf(fp, "%d", &rows);
		fgetc(fp);
		fgetc(fp);
		fscanf(fp, "%d", &cols);
		if (cols > length) cols = length;
		fgetc(fp);
		fgetc(fp);
		fscanf(fp, "%d", &udate);
		fgetc(fp);
		fgetc(fp);
		fgetc(fp);
		p = (node*)malloc(sizeof(node));
		p->col = cols;
		p->row = rows;
		p->date = udate;
		q->next = p;
		q = p;
	}
	q->next = NULL;
	q = head;
	head = head->next;
	free(q);
	return head;
}
void destroyList(node* head) 
{
	node* p, * q;
	for (p = head; p != NULL;)
	{
		q = p;
		p = p->next;
		free(q);
	}
}
char* my_itoa(char* str, int num)//局部变量，出了作用域释放
{
	int i = 0;
	while (num != 0)
	{
		str[i] = num % 10 + '0';//得到字符串“3 2 1”
		num = num / 10;//num == 0
		i++;//i == 3
	}
	str[i] = '\0';//str[3]
	i--;
	for (int j = 0; j < i; j++, i--)
	{
		char ch = str[j];//0号下标的元素赋给ch
		str[j] = str[i];//2号下标的元素赋给0号下标==》str【】==121
		str[i] = ch;//把3赋给2号下标
	}
	return str;
}
void suanfa(node* head1, node* head2)
{
	//FILE* fp;
	int i，j;
	float sum;
	node* p, * q;
	char ss[10];
		//使用遗传算法内容
	  {

		for (int i = 1;i <= length;i++) {
			X[i-1] = gen_group[0].info[i];
		}

	  }
	  printf("\nit is %f\n", gen_group[0].suitability);
		 char s[8] = "";

		 FILE* fp;//遗传算法
		 fp = fopen("X.txt", "w");
		 if (fp == NULL)
		 {
			 printf("文件打开失败\n");
		 }
		 else
		 {
			 for (int i = 0; i < length; i++)
			 {
				 //fputs(my_itoa(s, X[i]), fp);
				 sprintf(ss, "%f", X[i]); // double 到 char 
				 fputs(ss, fp);
				 fputs("\n", fp);
			 }
		 }
		 fclose(fp);

		 float sum2 = 0;
		 float  max = 0;
		 int j;
		 int  i = j = 1;
		 p = head1;
		 //printf("19");
		 q = head2;
		 while (q != NULL && i < length) {
			 while (p->row == i && p != NULL) {		//按行进行计算  第i列对应x的第p->col行    乘积相加  再将结果减去B对应的值   q->date  //i<=196608
				 sum2 = sum2 + 1.0 * p->date * 1.0 * X[p->col -1] * 1.0;
				 p = p->next;
			 }
			 i++;//下一行
			 sum2 -= (q->date) * 1.0;//AX-B b的值是依次排列的
			 if (sum2 < 0) {//绝对值
				 sum2 = -sum2;
				 printf(" \n sum2 %f \n", sum2);
			 }
			 q = q->next;
			 printf(" \n sum2 %f \n", sum2);
			 if (sum2 > max) {
				 max = sum2;//max的理想状态应该是0
				 //printf("\n\n max  %f   \n\n", max);
			 }
			 sum2 = 0; //重置sum2计算下一个
			 if (p == NULL || q == NULL) {
				 break;
			 }
		 }
		 printf("\n\n max  %f   \n\n", max);
		//不使用遗传算法的内容
		  FILE* fp2;
		  fp2 = fopen("X2.txt", "w");
		memset(X, 0, sizeof(X));
		memset(a, 0, sizeof(a));
		p = head1;
		float flag;
		q = head2;
		for (int j = 0; j <length;j++)
		{
			sum = 0;
			while (p->row == j + 1)
			{
				if (j + 1 == p->col)
				{
				flag =1.0* p->date;
					p = p->next;
				}
				else
				{
					sum += (-1) * p->date * X[p->col - 1];
					p = p->next;
				}
				if (p == NULL)
					break;
			}
			a[j] = X[j] = (sum + q->date) / flag;
			q = q->next;
		}
		if (fp2 == NULL)
		{
			printf("文件打开失败\n");
		}
		else
		{
			for (int i = 0; i < length; i++)
			{
				sprintf(ss, "%f", X[i]); // double 到 char 
				fputs(ss, fp2);
				fputs("\n", fp2);
			}
		}
		fclose(fp2);
		p = head1;
		q = head2;
		max = 0;
		i = j = 1;
		sum2 = 0;
		printf("不使用遗传算法");
		while (q != NULL && i < length) {
			while (p->row == i && p != NULL) {		//按行进行计算  第i列对应x的第p->col行    乘积相加  再将结果减去B对应的值   q->date  //i<=196608
				sum2 = sum2 + 1.0 * p->date * 1.0 * X[p->col - 1] * 1.0;
				p = p->next;
			}
			i++;//下一行
			sum2 -= (q->date) * 1.0;//AX-B b的值是依次排列的
			if (sum2 < 0) {//绝对值
				sum2 = -sum2;
				printf(" \n sum2 %f \n", sum2);
			}
			q = q->next;
			printf("\n  sum   %f  \n", sum2);
			if (sum2 > max) {
				max = sum2;//max的理想状态应该是0
				//printf("\n\n max  %f   \n\n", max);
			}
			sum2 = 0; //重置sum2计算下一个
			if (p == NULL || q == NULL) {
				break;
			}
		}
		printf("\n\n max  %f   \n\n", max);
}

使用该方法
达到的精确度可以到小数点后两位。