修正的高斯-牛顿法

如果已知探索序列的某个点Kn，在Kn附近把函数fl近似当作线性函数［3，31］，于是可得在Kn的泰勒展开式

含水层参数识别方法

式中ΔK=K-Kn，ΔKj（j=1，2，…，m）是它的分量。

将式（4-4）代入式（4-1），得

含水层参数识别方法

式中L是观测数据总数。

按照极值的必要条件，应有

含水层参数识别方法

将式（4-5）代入式（4-6），得

含水层参数识别方法

整理上式，得

含水层参数识别方法

以上方程组可表示为以下形式

含水层参数识别方法

式中矩阵

含水层参数识别方法

为L×M阶矩阵。

由式（4-8），得

含水层参数识别方法

于是可得改进点

含水层参数识别方法

将上式修正，得如下修正的高斯-牛顿法迭代序列

含水层参数识别方法

式中-Pn 和λn 分别是第n 次迭代的探索方向和步长。如果已知初始估计值 K0，则由上式可以算出 K1。重复这一过程，直到满足预先给定的迭代次数 N 或计算的Kn+1使得目标函数的绝对值小于ζ，即E（K）≤ζ，这里ζ是一个小的正数。为了确定方向 Pn ，首先需要计算目标函数 E 的梯度g 和高斯-牛顿方向ΔK。这样

含水层参数识别方法

式中：

含水层参数识别方法

式中[J]是灵敏度系数矩阵。现在引入投影算子-P 以保证Kn+1在给定的约束范围内，-P 的定义为

含水层参数识别方法

方程（4-10）中的步长是一个标量，可以用二次插值公式确定。

图4-1 用修正的高斯-牛顿法识别含水层参数的流程图

#include<stdio.h>

#include<math.h>

double f(double x)

{

    return (x*x*x-3*x-1)

}

double dfdx(double x)

{

    return (3*x*x-3)

}

void main()

{

    double x0=2,x1,y,err

    int n=0

    double tol=0.5*1.0e-3

    printf("解议程x*x*x-3*x-1=0的初始值为x0=%f,误差限为%f\n",x0,tol)

    do{

        y=dfdx(x0)

        if(fabs(y)<1.0e-5){

            printf("奇异:函数在%f的导数接近0",x0)

            break

        }else{

            x1=x0-f(x0)/y

            err=x0-x1

            n++

            printf("迭代%d次的解x%d=%f,error=%f,函数值f(x%d)=%f\n",n,n,x1,err,n,f(x1))

            x0=x1

        }

    }while(fabs(err)>tol)

}

牛顿-拉夫逊法早在50年代末就已应用于求解电力系统潮流问题，但作为一种实用的，有竞争力的电力系统潮流计算方法，则是在应用了稀疏矩阵技巧和高斯消去法求修正方程后。牛顿-拉夫逊法是求解非线性代数方程有效的迭代计算。

1.3MATLAB概述

目前电子计算机已广泛应用于电力系统的分析计算，潮流计算是其基本应用软件之一。现有很多潮流计算方法。对潮流计算方法有五方面的要求：（1）计算速度快（2）内存需要少（3）计算结果有良好的可靠性和可信性（4）适应性好，亦即能处理变压器变比调整、系统元件的不同描述和与其它程序配合的能力强（5）简单。

MATLAB是一种交互式、面向对象的程序设计语言，广泛应用于工业界与学术界，主要用于矩阵运算，同时在数值分析、自动控制模拟、数字信号处理、动态分析、绘图等方面也具有强大的功能。

MATLAB程序设计语言结构完整，且具有优良的移植性，它的基本数据元素是不需要定义的数组。它可以高效率地解决工业计算问题，特别是关于矩阵和矢量的计算。MATLAB与C语言和FORTRAN语言相比更容易被掌握。通过M语言，可以用类似数学公式的方式来编写算法，大大降低了程序所需的难度并节省了时间，从而可把主要的精力集中在算法的构思而不是编程上。

另外，MATLAB提供了一种特殊的工具：工具箱（TOOLBOXES）.这些工具箱主要包括：信号处理（SIGNAL PROCESSING）、控制系统（CONTROL SYSTEMS）、神经网络（NEURAL NETWORKS）、模糊逻辑(FUZZY LOGIC)、小波(WAVELETS)和模拟（SIMULATION）等等。不同领域、不同层次的用户通过相应工具的学习和应用，可以方便地进行计算、分析及设计工作。

MATLAB设计中，原始数据的填写格式是很关键的一个环节，它与程序使用的方便性和灵活性有着直接的关系。

原始数据输入格式的设计，主要应从使用的角度出发，原则是简单明了，便于修改。

2．1 电力系统的基本概念

2．1．1电力系统

（1）电力系统：发电机把机械能转化为电能，电能经变压器和电力线路输送并分配到用户，在那里经电动机、电炉和电灯等设备又将电能转化为机械能、热能和光能等。这些生产、变换、输送、分配、消费电能的发电机、变压器、变换器、电力线路及各种用电设备等联系在一起组成的统一整体称为电力系统。

（2）电力网：电力系统中除发电机和用电设备外的部分。

（3）动力系统：电力系统和“动力部分”的总和。

“动力部分”：包括火力发电厂的锅炉、汽轮机、热力网和用电设备，水力发电厂的水库和水轮机，核电厂的反应堆等。

2．1．2电力系统的负荷和负荷曲线

（1）电力系统的负荷：系统中千万个用电设备消费功率的总和，包括异步电动机、同步电动机、电热炉、整流设备、照明设备等若干类。

（2）电力系统的供电负荷：综合用电负荷加上电力网中损耗的功率。

（3）电力系统的发电负荷：供电负荷加上发电厂本身的消耗功率。

（4）各用电设备的有功功率和无功功率随受电电压和系统频率的变化而变化，其变化规律不尽相同，综合用电负荷随电压和频率的变化规律是各用电负荷变化规律的合成。

（5）负荷曲线：某一时间段内负荷随时间而变化的规律。

（6）按负荷种类可分有功功率负荷和无功功率负荷；按时间长短可分为日负荷和年负荷曲线；按计量地点可分为个别用户、电力线路、变电所、发电厂以至整个系统的负荷曲线。将上述三种分类相结合，就确定了某一种特定的负荷曲线。不同行业的有功功率日负荷曲线差别很大。负荷曲线对电力系统的运行又很重要的意义，它是安排日发电计划，确定各发电厂发电任务以及确定系统运行方式等的重要依据。

2．2 电力系统的基本元件

2．2．1 发电机

现代电力工业中，无论是火力发电、水力发电或核能发电，几乎全部采用同步交流发电机。电机的电枢布置在定子上，励磁绕组布置在转子上，作为旋转式磁极。同步发电机的转速（转/MIN）和系统频率f（HZ）之间有着严格的关系，即n=60f/p式中p为电机的极对数。

根据转子结构型式的不同，分为隐极式和凸极式发电机，前者转子没有显露出来的磁极，后者则有。

转子的励磁型式有直流励磁系统和可控硅励磁系统，后者利用同轴交流励磁机或由同步发电机本身发出的交流电，经整流后供给转子。直流励磁机有换向问题，故其制造容量受到限制，所以，在大容量发电机中均可采用可控硅励磁系统。

2．2．2 电力变压器

电力变压器是电力系统中广泛使用的升压和降压设备。据统计，电力系统中变压器的安装总容量约为发电机安装容量的6－8倍。按用途，电力变压器可分为升压变压器、降亚变压器、配电变压器和联络变压器。按相数分，变压器可分为单相式和三相式。按每相线圈分，又有双绕组和三绕组之分。按线圈耦合的方式，可分为普通变压器和自耦变压器。

2．2．3 电力线路

（1）架空线路：由导线、避雷针、杆塔、绝缘子和金具等构成。

（2）电缆线路：由导线、绝缘层、包护层等构成。

2．2．4 无功功率补偿设备

主要的无功功率补偿设备有同步调相机、电力电容器和静止补偿器。

2．3 电力系统元件的数学模型

2．3．1 电力线路的等值电路

在电力系统分析中，一般只考虑电力线路两侧端口的电压和电流，把电力线路作为无源双口网络处理。

线路的双口网络方程：

Z=B=*L*

2．3．2 变压器的等值电路

（1）双绕组变压器等值电路

（2）三绕组变压器等值电路

2．3．3 同步发电机的数学模型

2．3．4 电力系统负荷

2．3．5 多级电压电力系统的等值电路

2．4 电力系统稳态运行分析

2．4．1 电力线路的电压损耗与功率损耗

2．4．2 变压器中的功率损耗与电压损耗

2．4．3 辐射形网络的分析计算

辐射形电力网的特点是各条线路有明确的始端与末端。辐射形电力网的分析计算就是利用已知的负荷、节点电压来求取未知的节点电压、线路功率分布、功率损耗及始端输出功率。

辐射形电力网的分析计算，根据已知条件的不同分两种

1 已知末端功率与电压：即 从末端逐级往上推算，直至求得各要求的量。

2 已知末端功率、始端电压：末端可理解成一负荷点，始端为电源点或电压中枢点。采用迭代法。

（1）假设末端电压为线路额定电压，利用第一种方法求得始端功率及全网功率分布。

（2）用求得的线路始端功率和已知的线路始端电压，计算线路末端电压和全网功率分布。

（3）用第（2）步求得的线路末端电压计算线路始端功率和全网功率分布，如求得的各线路功率与前一次相同计算的结果相差小于允许值，就可以认为本步求得的线路电压和全网功率分布为最终计算结果。否则，返回第二步重新进行计算。

2．4．4 复杂电力系统潮流计算

电力系统潮流计算始对复杂电力系统正常和故障条件下稳态运行状态的计算。潮流计算的目标始求取电力系统在给定运行方式下的节点电压和功率分布，用以检查系统各元件是否过负荷、各点电压是否满足要求、功率的分布和分配是否合理以及功率损耗等。对现有电力系统的运行和扩建，对新的电力系统进行规划设计以及对电力系统进行静态和暂态稳定分析都是以潮流计算为基础。因此，潮流计算是电力系统计算分析中的一种最基本的计算。

潮流计算结果的用途，例如用于电力系统稳定研究、安全估计或最优潮流等也对潮流计算的模型和方法有直接影响。

2．5 电力系统潮流计算机算法

2．5．1电力系统潮流计算机算法概述

2．5．1．1 导纳矩阵的形成

2．5．1．2 节点类型

（1）PV节点：柱入有功功率P为给定值，电压也保持在给定数值。

（2）PQ节点：诸如有功功率和无功功率是给定的。

（3）平衡节点：用来平衡全电网的功率。选一容量足够大的发电机担任平衡全电网功率的职责。

平衡节点的电压大小与相位是给定的，通常以它的相角为参考量，即取其电压相角为0。一个独立的电力网中只设一个平衡点。

2．5．1．3 高斯迭代法

2．5．2 牛顿－拉夫逊法

2．5．2．1 原理

2．5．2．2 基本步骤

基本步骤：

（1）形成节点导纳矩阵

（2）将各节点电压设初值U，

（3）将节点初值代入相关求式，求出修正方程式的常数项向量

（4）将节点电压初值代入求式，求出雅可比矩阵元素

（5）求解修正方程，求修正向量

（6）求取节点电压的新值

（7）检查是否收敛，如不收敛，则以各节点电压的新值作为初值自第3步重新开始进行狭义次迭代，否则转入下一步

（8）计算支路功率分布，PV节点无功功率和平衡节点柱入功率。

2．5．2．3 注意事项

2．5．2．4 程序流程框图

2．6 软件设计

2．6．1 方案选择及说明

2．6．2 方案求解

2．6．3 MATLAB编程说明及元件描述

2．6．4 程序

#include<stdio.h>

struct powernode

{

float pi

float qi

int i

float vi

}

struct powernode wg[20]

struct powernode wl[20]

struct linedata

{

int i

int j

float r

float x

float y /*包括变压器变比*/

float k /*只用作标析变压器，变压器变比仍在y中*/

}

struct linedata zl[20]

struct linedata t3/*临时数组*/

static double y[][3]/*在matrixform中应用*/

int t=0

int t2,ti,tj /*临时记数单元*/

float temp

float tx,tr,YK /*中间工作单元(在matrixform中应用)*/

double GIJ,BIJ/*中间工作单元(在matrixform中应用)*/

int N /*总节点数*/

int zls

int Q,V,PVS,PVD

int GS

int LS

float vo

float Eps

static double GII[]={0},BII[]={0},YDS[]={0},YDZ[]={0},B[]={0}/*添加数组*/

/*因子表形成时定义的数据*/

struct pvdata

{

float vis

int i

}

static struct pvdata pv[]={0}

datain()

{

clrscr()

printf("program runningn" )

printf("n")

printf("please input the aggregate to the system note")/*总节点数*/

scanf("%d",&N)

printf("n")

printf(" PQ note IN ALL?")/*总节点数*/

scanf("%d",&Q)

PVS=(N-Q)-1

printf("n")

printf("them input the aggregate to the system power line")

scanf("%d",&zls)/*输电线路数和变压器的总数*/

printf("n")

printf("electromotor node in all :?")/*发电机节点总数*/

scanf("%d",&GS)

printf("n")

printf("load node in all : ?")/*负荷节点总数*/

scanf("%d",&LS)

printf("n")

printf("average electric voltage")/*平均电压*/

scanf("%f",vo)

printf("n")

printf("n")

printf("please input the date messagen")

printf("follow the format like it: i,j,r,x,y,kn")

do{

t++

scanf("%d,%d,%f,%f,%f",&zl[t].i,&zl[t].j,&zl[t].r,&zl[t].x,&zl[t].y,&zl[t].k)

printf("processing....n")

if(zl[t].i>zl[t].j)

{

temp=zl[t].i

zl[t].i=zl[t].j

zl[t].j=temp

/* if(zl[t].k!=1) */ /*要考虑归算问题不？？？？*/

}

printf("data you input is:n " )

printf("%d,%d,%f,%f,%f",zl[t].i,zl[t].j,zl[t].r,zl[t].x,zl[t].y,zl[t].k)

}while(zl[t].i!=0&&zl[t].j==0)

for(t2=tt>0t--) /*冒泡法排序*/

{

for(t2>0t2--)

{

if(zl[t2].i<zl[t2-1].i)

{

t3.i=zl[t2].it3.j=zl[t2].jt3.r=zl[t2].rt3.x=zl[t2].xt3.y=zl[t2].yt3.k=zl[t2].k

zl[t2].i=zl[t2-1].izl[t2].j=zl[t2-1].jzl[t2].r=zl[t2-1].rzl[t2].x=zl[t2-1].xzl[t2].y=zl[t2-1].yzl[t2].k=zl[t2-1].k

zl[t2-1].i=t3.izl[t2-1].j=t3.jzl[t2-1].r=t3.rzl[t2-1].x=t3.xzl[t2-1].y=t3.yzl[t2-1].k=t3.k

}

else if(zl[t2].i==zl[t2-1].i)

{if(zl[t2].j<zl[t2-1].j)

{

t3.i=zl[t2].it3.j=zl[t2].jt3.r=zl[t2].rt3.x=zl[t2].xt3.y=zl[t2].yt3.k=zl[t2].k

zl[t2].i=zl[t2-1].izl[t2].j=zl[t2-1].jzl[t2].r=zl[t2-1].rzl[t2].x=zl[t2-1].xzl[t2].y=zl[t2-1].yzl[t2].k=zl[t2-1].k

zl[t2-1].i=t3.izl[t2-1].j=t3.jzl[t2-1].r=t3.rzl[t2-1].x=t3.xzl[t2-1].y=t3.yzl[t2-1].k=t3.k

}

}

}

}

printf("n")

t=0

printf("please input wg~!n")

do

{

scanf("%f,%f,%d,%f",&wg[t].pi,&wg[t].qi,&wg[t].i,&wg[t].vi)

t++

}while(t!=GS)ti=0

for(t=0t<GSt++){if(wg[t].vi<o){pv[ti].vis=labs(wg[t].vi)pv[ti].i=wg[t].iti++}}

t2=0

printf("please input WL~!n")

do

{

scanf("%f,%f,%d,%f",&wl[t2].pi,&wl[t2].qi,&wl[t2].i,&wl[t2].vi)

t2++

}while(t2!=LS)

for(t=0t<LSt++){if(wl[t].vi<o){pv[ti].vis=labs(wl[t].vi)pv[ti].i=wl[t].iti++}}

}

matrixform()

{

for(t=1t<Nt++)

{

GII[t]=0

BII[t]=0

YDS[t]=0

}

for(t2=1t<zlst2++)

{

ti=labs(zl[t2].i)

tj=labs(zl[t2].j)

tr=zl[t2].r

tx=zl[t2].x

temp=ldexp(tr,1)+ldexp(tx,1)

GIJ=tr/tempBIJ=tx/temp

y[t2][1]=-GIJ

y[t2][2]=-BIJ

y[t2][3]=tj

GII[ti]=GII[ti]+GIJBII[ti]=BII[ti]+BIJ

GII[tj]=GII[tj]+GIJBII[tj]=BII[tj]+BIJ

YDS[ti]=YDS[ti]+1

}

YDZ[1]=1

for(t=1t<N-1t++)

{

YDZ[t+1]=YDZ[t]+YDS[t]

} /*矩阵型成第一部完成*/

/*矩阵型成第二部开始*/

for(t2=1t<zlst2++)

{ /*.k只用作变压器的标析，变压器变比仍在y中*/

ti=zl[t2].itj=zl[t2].jYK=zl[t2].y

if(ti<0||tj<0)

{ if(ti<0)

ti=labs(ti)

else

ti=labs(tj)

GIJ=y[t2][1]BIJ=y[t2][2]

GII[t2]=GII[t2]+(1-1/YK/YK)*GIJ

BII[t2]=BII[t2]+(1-1/YK/YK)*BIJ

y[t2][1]=GIJ/YK

y[t2][2]=BIJ/YK

}

else

GIJ=0

BIJ=YK/2

SY(tr) /*这个东东要调用，实现节点累计自导纳*/

SY(tj) /*SY的过程是完成向一个节点累计相应自导纳的实部和虚部*/

}

}

int sign,ld,k2,x,im,ai /*k2控制台开关，负荷静态特性开关*/

static float fd[]={0}

unsigned AF[1]

static int u[]={0}/*???????????怎么实现？来自那里???????*/

divisorform()

{

/*暂时不知道LD PVD 等的作用……待善*/

PVD=pv[0].i

ld=wl[0].i

t=0

do{

t2++

if(sign==1&&t2==PVD)

{t=t+1pvd=pv[t].ifd[t2]=0di[t2]=0

if(k2==0&&t2==ld)

{t2=t2+1ld=wl[t2].i}

}continue

else

B[t2]=BII[t2]

if(k2==0&&sing==1&&t2==ld)

{

B[t2]=B[t2]+AF[1]*wl[t2].qi/wl[t2].vi/wl[t2].vit2=t2+1ld=wl[t2].i

}

for(temp=YDZ[t2]temp<YDZ[t2+1]-1temp++)

{

tj=Y[temp][3]B[tj]=Y[temp][2]

}

if(sign=1)

{for(temp=1temp<PVStemp++)

tj=pv[temp][2]

B[tj]=0

}

x=2im=1

do{im++

if(im>t2-1)

break

else

temp=1

for(temp!>fd[im]){if(u[x+1]!=1){temp=temp+1x=x+2}else ai=u[x]/} /*u[]未完成*/

continue

}

}while(t2!=N-1)

}

dataout()

{

clrscr()

printf("note 1 voltagen")

printf("(.639696730300784) + j (1.832939) = 1.94136001255537 ∠ 70.7609880529659°n")

printf("87u&婾[1]??u?孢???�u

--------------------------------------------------------------------------------

??虍鉧C&8u謤蛝髻??n")

}

main()

{

datain() /*数据输入及处理*/

matrixform() /*矩阵的形成*/

/* divisorform() */ /*因子表的形成*/

matrixsolve() /*矩阵线形方程的求解*/

/* nodepower()*/ /*迭代过程中节点功率的计算*/

/* iterate() */ /*迭代*/

dataout() /*数据输出及支路功率计算*/

}

---