基于C++的PQ分解法求解潮流计算1——高斯消元法

努力working的kk

已于 2023-05-01 14:14:22 修改

阅读量1.4k

点赞数 3

分类专栏：电力系统分析文章标签： c++ 经验分享课程设计 gitcode

于 2023-04-28 14:30:56 首次发布

本文链接：https://blog.csdn.net/q158369751/article/details/130395686

版权

电力系统分析专栏收录该内容

2 篇文章 2 订阅

订阅专栏

系列文章目录

第一章使用c++实现PQ分解法1——高斯消元法
第二章基于C++的PQ分解法求解潮流计算2——因子表法

文章目录

系列文章目录
前言
一、高斯消元法简介
二、关键函数
三、总结

前言

从本章开始，我们将逐步讲解如何使用c++语言实现PQ分解法，并将其应用到电力系统的潮流计算之中。虽然是基于C++语言的代码，但绝大部分用的都是c语言的特性。
PQ分解法是依据高电压等级电力系统特性对牛顿拉夫逊方法进行改造而获得的方法，是用来求解潮流方程（多元非线性方程组）的经典算法，本系列文章将会对其的c++实现方法进行详细介绍。
具体而言，本系列文章需要实现以下任务。
（1）应用高斯消元法求解一次多元线性方程组
（2）应用因子表法完成重复求解同一系数矩阵的不同多元线性方程组
（3）使用稀疏技术和稀疏相量法对求解多元线性方程组的过程进行简化
（4）使用PQ分解法将多元非线性方程组转化为多元线性方程组
（5）使用节点编号优化法对PQ分解法进行简化
（6）将上述五个代码合并为一个文件，使用算例对其效果进行验证
本章则将介绍应用高斯消元法求解线性方程组的全部代码过程。另外笔者是第一次写文章，而且是c++的初学者，在实现代码过程中，遇到了很多非常弱智的语法错误和小坑，都记录下来留给自己借鉴。望不要苛责。

一、高斯消元法简介

高斯消元法是我们初中时就已经学习过的求解一次多元线性方程组的方法，思路非常简单，对于小规模问题很容易手算出来。对于中等规模的问题，可以使用线性代数的手段直接对增广矩阵进行处理，也可手算出来。但在电力系统中，动辄成千上万个节点，手算如此规模的线性方程组显然不大可能，所以需要使用计算机进行计算。高斯消元法虽然原理简单，但在代码实现过程中可能会遇到诸多问题，所以本文将尽可能详细地逐步讲解高斯消元法的所有细节，供读者参考。

二、关键函数

本节将对实现高斯消元法的所有函数进行讲解。若要使用c++实现高斯消元法，有如下几个函数是不可或缺的。
（1）互换一维数组
本函数可将两个一维数组的元素进行互换。需要注意边界条件，若输入的维数小于1则应报错。
（2）打印数组内容
为了方便调试，需要编写打印数组的函数。由于本实验中数组包含一维和二维数组，所以应当使用c++函数的多态特性，用一个函数满足打印一维和二维数组的需求。同时要注意在输出的合适位置添加空格和换行符，以提升可读性
（3）将系数矩阵转化为上三角矩阵
本函数为本文难点。首先要考虑每一行的第一个元素是否为零，如果为零则需要将非零的行替换上来，而且每迭代一次都需要考虑这个问题。其次是本文的各矩阵都是double类型，所以需要特殊的手段来判断一个值是否为零。最后就是在写函数时，要找清楚每个for循环的起始和终止条件，这个比较困难。另外需要满足可以处理任意维度二维数组的需求。
（4）根据上三角矩阵求解方程组
本函数同样为本文难点。判断每个for循环的起始和终止条件也需要反复琢磨。

1.互换一维数组

简单说一下笔者在这一函数所遇到的小坑。
①笔者一开始试图通过互换两个一维数组的首地址，从而完成数组元素互换的目的。这一设想对普通的数组适用，但用于互换二维数组内的两个行向量时出现了问题，并没有互换成功。原因是二维数组的每个元素都是按顺序在内存中排列的，无法直接将两个元素或行向量的地址进行互换或修改。
②然后笔者试图使用中间数组的方式对两个数组进行互换，但由于数组的维数不定，使用函数传参的话，无法在函数内创建中间数组，因为创建数组时必须使用常量而非变量。所以本法也不行。
③最后只能使用最简单的方法，通过循环的方式一个一个地用中间变量互换元素。

代码如下：

/*
* 本函数可将两个一维数组进行互换
* 本函数包含两个参数，分别是两个需要互换的数组名
*/
void swap_arrays(double* array1, double* array2, int n)
{
	
	double temp;
	int i = 0;
	if (n <= 0)
	{
		printf("the input array is wrong!");
	}
	else
	{
		for (i = 0; i < n; i++)
		{
			temp = array1[i];
			array1[i] = array2[i];
			array2[i] = temp;
		}
	}

}

2.打印数组内容

打印数组同样有很多小坑，主要是笔者写惯了c，很多c中的特性在c++中无法使用，所以出现了诸多问题。本函数最大的问题在于二维数组作为参数传递给函数时遇到的困难。因为在c++或c中，二维数组作为参数必须给出列数，但很多时候我们并不知道数组的列数，所以只能使用强制类型转换将二维数组的首地址传入函数。这样就使得我们必须在主调函数对二维数组名进行强制类型转换，这样就对使用该函数的用户带来了麻烦。但如果不使用强制类型转换，笔者也实在不知道如何操作。

代码如下：

/*
* 本函数可打印1、2维数组
* 本函数使用了多态特性，包括二或三个参数。打印一维数组的参数是数组名、列数，打印二维数组的参数是数组名、行数、列数，
* 二维数组作为参数传入本函数时，需要使用强制类型转换，将二维数组名转换为普通的地址。如print_array((double*)test1,2,3);
*/
void print_array(double *array, int a, int b)
{
		for (int i = 0; i < a; i++)
		{
			for (int j = 0; j < b; j++)
			{
				printf("%lf ", *(array + b * i + j));
			}
			printf("\n");
		}
}

void print_array(double* array,  int a)
{

		for (int i = 0; i < a; i++)
		{
			printf("%lf ", *(array + i));
		}
		printf("\n");
}

3.将系数矩阵转化为上三角矩阵

①本函数为实现能处理任意行列数的二维矩阵，采用了动态二维数组的处理方式
②由于二维数组中的元素均为double类型，所以定义zero = 1e-6，并认为所有小于zero的值均为0
③在将系数矩阵转化为上三角矩阵时，需迭代使用初等行变换。并要求在每次迭代时，不能出现某一列对角元素为0，而其下面的元素存在非0的情况。所以函数中第一个双for循环便是为了处理这种问题。
④函数中第二个双for循环是逐行进行初等行变换，消去列元素的操作。需要重点关注初始条件的设置方法

代码如下：

/*
* 本函数可将增广矩阵中的系数矩阵转化为上三角矩阵，同时初等行变换也会作用于常数列，最终获得的二维数组是原增广矩阵同解方程组的增广矩阵
* 本函数包括3个参数，分别是二维数组（增广矩阵）的动态数组形式，即指向指针的指针，指针中每一个元素都对应二维数组的一行的首地址。第2和3参数是二维数组的行和列数
* 在使用本函数时，需在主调函数中将二维数组名转化为指针形式，也就是将二维数组首地址传入函数，这样是为了解决将任意行列数的二维数组传入本函数的函数泛用性问题
*/
void matrix2triangular(double* arr, int n, int b)
{
		int i = 0;  // i为迭代次数
	int j = 0;  // j为迭代过程中临时变量
	int k = 0;  // k为迭代过程中临时变量
	double factor = 0;  // factor为初等行变换时所乘的因数
	
	// 使用二维动态数组方式处理增广矩阵，将形参arr内容全部赋值给aug_matrix，对aug_matrix的处理同时会改变arr，而且将前者释放并不会影响arr
	double** aug_matrix = new double* [n];
	for (int i = 0; i < n; i++) {
		aug_matrix[i] = arr+b*i;
	}

	// 将本次迭代需要处理的所有行向量中首元素为非零的置于第一行
	for (i = 0; i < n; i++)
	{
		for (j = i + 1; j < n; j++)
		{

			if (fabs(aug_matrix[i][i]) < ZERO)
			{
				swap_arrays(aug_matrix[i], aug_matrix[j], n + 1);
			}
			else { break; }
		}

		for (j = i + 1; j < n; j++)
		{
			if (fabs(aug_matrix[j][i]) < ZERO)
				continue;

			else
			{
				factor = aug_matrix[j][i] / aug_matrix[i][i];  //factor中必须是aug_matrix[j][i]，不能是aug_matrix[i + 1][i]
				for (k = i; k < n + 1; k++)
				{
					aug_matrix[j][k] -= factor * aug_matrix[i][k];
				}
			}
		}
	}
	for (int i = 0; i < 3; i++) {
		for (int j = 0; j < 4; j++) {
			cout << aug_matrix[i][j] << " ";
		}
		cout << endl;
	}
	refine_arrays((double*)aug_matrix[0], n, n + 1);  // 将二维数组中过小的元素置0
	delete[] aug_matrix;
}

4.根据上三角矩阵求解方程组

本函数可以求取系数矩阵为上三角矩阵的增广矩阵的解。如double aug_matrix[3][4] = {1,1,1,12,1,1,0,10,1,2,5,22}；解得三个变量的值分别为8、2、2.
本函数用法与matrix2triangular完全一致

代码如下：

/*
* 本函数可对系数矩阵为上三角矩阵的增广矩阵进行求解，得到常数列，也就是每个方程组每个x的解
* 本函数包括3个参数，分别是二维数组（增广矩阵）的动态数组形式，即指向指针的指针，指针中每一个元素都对应二维数组的一行的首地址。第2和3参数是二维数组的行和列数
* 在使用本函数时，需在主调函数中将二维数组名转化为指针形式，也就是将二维数组首地址传入函数，这样是为了解决将任意行列数的二维数组传入本函数的函数泛用性问题
*/
void triangular_soving(double* arr, int n, int b)
{
	int i = 0;
	int j = 0;
	// 使用二维动态数组方式处理增广矩阵，将形参arr内容全部赋值给aug_matrix，对aug_matrix的处理同时会改变arr，而且将前者释放并不会影响arr
	double** aug_matrix = new double* [n];
	for (i = 0; i < n; i++) {
		aug_matrix[i] = arr + b * i;
	}
	//通过行循环和列循环对常数列逐个运算
	for (i = n - 1; i >= 0; i--)
	{
		for (j = n - 1; j >= i + 1; j--)
		{
			if (fabs(aug_matrix[j][j]) < ZERO) { continue; }
			else {
				aug_matrix[i][n] -= aug_matrix[i][j] / aug_matrix[j][j] * aug_matrix[j][n];
			}
		}
		aug_matrix[i][n] /= aug_matrix[i][i];
		aug_matrix[i][i] /= aug_matrix[i][i];
	}
	refine_arrays((double*)aug_matrix[0], n, n + 1);  // 将二维数组中过小的元素置0
	delete[] aug_matrix;
}

5.全部代码

下面为高斯消元法求解多元线性方程组的全部代码，可以用c++直接运行

#include <stdio.h>
#include <math.h>
#include <iostream>
using namespace std;
constexpr auto ZERO = 1e-6;;

/*
* 本函数可将两个一维数组进行互换
* 本函数包含两个参数，分别是两个需要互换的数组名
*/
void swap_arrays(double* array1, double* array2, int n)
{
	
	double temp;
	int i = 0;
	if (n <= 0)
	{
		printf("the input array is wrong!");
	}
	else
	{
		for (i = 0; i < n; i++)
		{
			temp = array1[i];
			array1[i] = array2[i];
			array2[i] = temp;
		}
	}

}

/*
* 本函数可打印1、2维数组
* 本函数使用了多态特性，包括二或三个参数。打印一维数组的参数是数组名、列数，打印二维数组的参数是数组名、行数、列数，
* 二维数组作为参数传入本函数时，需要使用强制类型转换，将二维数组名转换为普通的地址。如print_array((double*)test1,2,3);
*/
void print_array(double *array, int a, int b)
{
		for (int i = 0; i < a; i++)
		{
			for (int j = 0; j < b; j++)
			{
				printf("%lf ", *(array + b * i + j));
			}
			printf("\n");
		}
}

void print_array(double* array,  int a)
{

		for (int i = 0; i < a; i++)
		{
			printf("%lf ", *(array + i));
		}
		printf("\n");
}

/*
* 本函数可将二维数组小于ZERO的值置零
* 本函数包含三个参数，二维数组的首地址和行数列数
*/
void refine_arrays(double* array, int a, int b)
{
	for (int i = 0; i < a; i++)
	{
		for (int j = 0; j < b; j++)
		{
			if (fabs(*(array + i * b + j)) < ZERO)
			{
				*(array + i * b + j) = 0.0;
			}
		}
	}
}

/*
* 本函数可将增广矩阵中的系数矩阵转化为上三角矩阵，同时初等行变换也会作用于常数列，最终获得的二维数组是原增广矩阵同解方程组的增广矩阵
* 本函数包括3个参数，分别是二维数组（增广矩阵）的动态数组形式，即指向指针的指针，指针中每一个元素都对应二维数组的一行的首地址。第2和3参数是二维数组的行和列数
* 在使用本函数时，需在主调函数中将二维数组名转化为指针形式，也就是将二维数组首地址传入函数，这样是为了解决将任意行列数的二维数组传入本函数的函数泛用性问题
*/
void matrix2triangular(double* arr, int n, int b)
{
	int i = 0;  // i为迭代次数
	int j = 0;  // j为迭代过程中临时变量
	int k = 0;  // k为迭代过程中临时变量
	double factor = 0;  // factor为初等行变换时所乘的因数
	
	// 使用二维动态数组方式处理增广矩阵，将形参arr内容全部赋值给aug_matrix，对aug_matrix的处理同时会改变arr，而且将前者释放并不会影响arr
	double** aug_matrix = new double* [n];
	for (int i = 0; i < n; i++) {
		aug_matrix[i] = arr+b*i;
	}

	// 将本次迭代需要处理的所有行向量中首元素为非零的置于第一行
	for (i = 0; i < n; i++)
	{
		for (j = i + 1; j < n; j++)
		{

			if (fabs(aug_matrix[i][i]) < ZERO)
			{
				swap_arrays(aug_matrix[i], aug_matrix[j], n + 1);
			}
			else { break; }
		}

		for (j = i + 1; j < n; j++)
		{
			if (fabs(aug_matrix[j][i]) < ZERO)
				continue;

			else
			{
				factor = aug_matrix[j][i] / aug_matrix[i][i];  //factor中必须是aug_matrix[j][i]，不能是aug_matrix[i + 1][i]
				for (k = i; k < n + 1; k++)
				{
					aug_matrix[j][k] -= factor * aug_matrix[i][k];
				}
			}
		}
	}
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n+1; j++) {
			cout << aug_matrix[i][j] << " ";
		}
		cout << endl;
	}
	refine_arrays((double*)aug_matrix[0], n, n + 1);  // 将二维数组中过小的元素置0
	delete[] aug_matrix;
}

/*
* 本函数可对系数矩阵为上三角矩阵的增广矩阵进行求解，得到常数列，也就是每个方程组每个x的解
* 本函数包括3个参数，分别是二维数组（增广矩阵）的动态数组形式，即指向指针的指针，指针中每一个元素都对应二维数组的一行的首地址。第2和3参数是二维数组的行和列数
* 在使用本函数时，需在主调函数中将二维数组名转化为指针形式，也就是将二维数组首地址传入函数，这样是为了解决将任意行列数的二维数组传入本函数的函数泛用性问题
*/
void triangular_soving(double* arr, int n, int b)
{
	int i = 0;
	int j = 0;
	// 使用二维动态数组方式处理增广矩阵，将形参arr内容全部赋值给aug_matrix，对aug_matrix的处理同时会改变arr，而且将前者释放并不会影响arr
	double** aug_matrix = new double* [n];
	for (i = 0; i < n; i++) {
		aug_matrix[i] = arr + b * i;
	}
	//通过行循环和列循环对常数列逐个运算
	for (i = n - 1; i >= 0; i--)
	{
		for (j = n - 1; j >= i + 1; j--)
		{
			if (fabs(aug_matrix[j][j]) < ZERO) { continue; }
			else {
				aug_matrix[i][n] -= aug_matrix[i][j] / aug_matrix[j][j] * aug_matrix[j][n];
			}
		}
		aug_matrix[i][n] /= aug_matrix[i][i];
		aug_matrix[i][i] /= aug_matrix[i][i];
	}
	refine_arrays((double*)aug_matrix[0], n, n + 1);  // 将二维数组中过小的元素置0
	delete[] aug_matrix;
}

int main()
{
	// 表示增广矩阵的行数
	int n = 0;
	int i = 0;
	int j = 0;
	double temp[100] = { 0 };
	cout << "please input the row of your augmented matrix" << endl;
	cin >> n;
	cout << "please input your augmented matrix" << endl;
	for (i = 0; i < n * (n + 1); i++)
	{
		cin >> temp[i];
	}
	cout << "your augmented matrix is as follows" << endl;
	print_array(temp, n, n + 1);
	cout << "your upper triangular augmented matrix is as follows" << endl;
	// 接下来的两个关键函数中的第一个参数都是二维数组的第一个元素的地址，所以如果用二维数组名传参，用（double*)强转就行。如果用一位数组名，则直接传递即可
	matrix2triangular(temp, n, n + 1);
	triangular_soving(temp, n, n + 1);
	cout << "your solution of augmented matrix is as follows" << endl;
	// 打印出方程组的解
	for (i = 1; i < n+1; i++)
	{
		cout << "x" << i <<"="<<temp[(n + 1) * i-1]<<endl;
	}
	
	return 0;
}

三、总结

以上便是基于c++的应用高斯消元法求解线性方程组的全过程。为了增加各函数的泛用性，很多关键函数关于二维数组传递的参数都采用了二维数组首行首元素地址的形式，所以只需要将增广矩阵按照顺序输入进一维数组，便可直接将一维数组名作为参数传给函数。因为一维数组名本身就是首元素的首地址。
本系列全部代码将会放入github中开源。

github链接: link