using System;
using System.Collections.Generic;
using System.Text;
using
using System.Collections;
using System.Timers;

namespace EletricFlow
{
    class Bus_type
    {
        public int type;
        public double vol, ang;
        public double p, q, q_vol;
        public double g0, b0;
        public int bus_serino;
        public Bus_type()
        {
            this.type = 0;
            this.vol = 0;
            this.ang = 0;
            this.p = 0;
            this.q = 0;
            this.q_vol = 0;
            this.g0 = 0;
            this.b0 = 0;
            this.bus_serino = 0;
        }
        public void clearall()
        {
            this.type = 0;
            this.vol = 0;
            this.ang = 0;
            this.p = 0;
            this.q = 0;
            this.q_vol = 0;
            this.g0 = 0;
            this.b0 = 0;
            this.bus_serino = 0;
        }
    }
    class Branch_type
    {
        public int from, to, start, end;
        public double r, x, b, k;
        public Complex Sij, Sji, delta_S;
        public Branch_type()
        {
            this.from = 0;
            this.to = 0;
            this.start = 0;
            this.end = 0;
            this.r = 0;
            this.x = 0;
            this.b = 0;
            this.k = 0;
            this.Sij = new Complex();
            this.Sji = new Complex();
            this.delta_S = new Complex();
        }
    }
    class GVariant
    {

        public ArrayList Bus = new ArrayList();////
        public ArrayList Branch = new ArrayList();
        public ArrayList NAM = new ArrayList();
        public int i, bus_no, branch_no, pq_no, pv_no, iter;
        public ArrayList Node_new = new ArrayList();
        public int[] Node_old;
        public int max_iter = 50;
        public GVariant()
        {
            i = 0;
            bus_no = 0;
            branch_no = 0;
            pq_no = 0;
            pv_no = 0;
            iter = 0;
        }


    }
    class FlowFunction
    {
        public static GVariant gv = new GVariant();
        public static Complex zero = new Complex(0, 0);/*常用的0和1*/
        public static Complex one = new Complex(1, 0);/*常用的0和1*/
        /*输入函数，从IEEE标准测试系统读入数据并赋值给变量*/
        public static void Input(string filename)/////
        {
            //定义变量
            char ctype = 'w';
            double sbase, p_load, q_load, p_gen, q_gen;

         

            //打开文件
            if (!File.Exists(filename))
            {
                Console.WriteLine("{0} does not exist.", filename);
                return;
            }
            else
            {
                StreamReader sr = new StreamReader(filename);
                string line;
                ArrayList tem_list = new ArrayList();
                line = sr.ReadLine();
                //Console.WriteLine(line);
                string[] line_a = line.Split(' ');
                for (int u = 0; u < line_a.Length; u++)
                {
                    if (line_a[u].Length != 0)
                        tem_list.Add(line_a[u]);
                }
                // Console.WriteLine(tem_list.Count);
                sbase = double.Parse((string)tem_list[1]);///读基准容量
                tem_list.Clear();

                line_a.Initialize();
                line = sr.ReadLine();

                line_a = line.Split(' ');
                for (int u = 0; u < line_a.Length; u++)
                {
                    if (line_a[u].Length != 0)
                        tem_list.Add(line_a[u]);
                }
                gv.bus_no = int.Parse((string)tem_list[3]);//母线节点数
                tem_list.Clear();
                line_a.Initialize();


                Bus_type[] Bus = new Bus_type[gv.bus_no + 1];
                for (int i = 0; i <= gv.bus_no; i++)
                    Bus[i] = new Bus_type();

                gv.Node_new.Add(0);
                gv.Node_old = new int[gv.bus_no + 1];
                for (int jj = 0; jj <= gv.bus_no; jj++)
                {
                    gv.Node_old[jj] = 0;
                }

                //   StreamWriter sw = new StreamWriter("new ieee300.txt");
                string[] bus_line = new string[gv.bus_no + 1];
                int max = 0;
                for (int i = 0; i <= gv.bus_no; i++)
                {
                    bus_line[i] = "";
                }
                for (int i = 1; i <= gv.bus_no; i++)
                {
                    line = sr.ReadLine();
                    bus_line[i] = line;
                    line_a = line.Split(' ');
                    for (int u = 0; u < line_a.Length; u++)
                    {
                        if (line_a[u].Length != 0)
                        {
                            tem_list.Add(line_a[u]);
                        }
                    }
                    if (tem_list.Count > max)
                        max = tem_list.Count;
                    tem_list.Clear();
                    line_a.Initialize();
                }

                //@@@@

                for (int i = 1; i <= gv.bus_no; i++)
                {
                    Bus[i].ang = 0;
                    line = bus_line[i];
                    if (line != null)
                    {
                        line_a = line.Split(' ');
                        for (int u = 0; u < line_a.Length; u++)
                        {
                            if (line_a[u].Length != 0)
                                tem_list.Add(line_a[u]);
                        }

                    }

                    //string sq = "";
                    //for (int u = 0; u < tem_list.Count; u++)
                    //{
                    //    //Console.Write(tem_list[u] + "##");
                    //    sq += tem_list[u] + "##";


                    //}
                    //sw.WriteLine(sq);
                    // Console.WriteLine("max:{0};tem_lis_count:{1}",max,tem_list.Count);
                    if (tem_list.Count < max)
                    {
                        ctype = '1';
                    }
                    else
                    {

                        if (((string)tem_list[5]).Length == 1)
                            ctype = char.Parse((string)tem_list[5]);

                    }



                    switch (ctype) //根据节点类型不同选择读入不同的信息
                    {

                        case '1': //PQ节点
                            gv.pq_no++;
                            Bus[gv.pq_no].type = 1;
                            Bus[gv.pq_no].vol = 1;
                            p_load = double.Parse((string)(tem_list[7]));
                            q_load = double.Parse((string)(tem_list[8]));
                            p_gen = double.Parse((string)(tem_list[9]));
                            q_gen = double.Parse((string)(tem_list[10]));
                            Bus[gv.pq_no].bus_serino = i;///add myself
                            Bus[gv.pq_no].p = (p_gen - p_load) / sbase;
                            Bus[gv.pq_no].q = (q_gen - q_load) / sbase;
                            Bus[gv.pq_no].g0 = double.Parse((string)tem_list[15]);
                            Bus[gv.pq_no].b0 = double.Parse((string)tem_list[16]);
                            gv.Node_new.Add(gv.pq_no);
                            gv.Node_old[gv.pq_no] = i;
                            break;
                        case '0':////PQ节点
                            gv.pq_no++;
                            Bus[gv.pq_no].type = 0;
                            Bus[gv.pq_no].vol = 1;
                            p_load = double.Parse((string)(tem_list[8]));
                            q_load = double.Parse((string)(tem_list[9]));
                            p_gen = double.Parse((string)(tem_list[10]));
                            q_gen = double.Parse((string)(tem_list[11]));
                            Bus[gv.pq_no].bus_serino = i;///add myself
                            Bus[gv.pq_no].p = (p_gen - p_load) / sbase;
                            Bus[gv.pq_no].q = (q_gen - q_load) / sbase;
                            Bus[gv.pq_no].g0 = double.Parse((string)tem_list[16]);
                            Bus[gv.pq_no].b0 = double.Parse((string)tem_list[17]);
                            gv.Node_new.Add(gv.pq_no);
                            gv.Node_old[gv.pq_no] = i;
                            break;
                        case '2': //PV节点
                            gv.pv_no++;
                            Bus[gv.bus_no - gv.pv_no].type = 2;
                            Bus[gv.bus_no - gv.pv_no].vol = double.Parse(((string)tem_list[6]));
                            p_load = double.Parse((string)(tem_list[8]));
                            p_gen = double.Parse((string)(tem_list[10]));
                            // Console.WriteLine("rrrrrr");
                            Bus[gv.bus_no - gv.pv_no].bus_serino = i;///add myself
                            Bus[gv.bus_no - gv.pv_no].p = (p_gen - p_load) / sbase;
                            Bus[gv.bus_no - gv.pv_no].q = 0;
                            Bus[gv.bus_no - gv.pv_no].vol = double.Parse(((string)tem_list[6]));
                            Bus[gv.bus_no - gv.pv_no].g0 = double.Parse((string)tem_list[16]);
                            Bus[gv.bus_no - gv.pv_no].b0 = double.Parse((string)tem_list[17]);
                            // Console.WriteLine("pvg0:{0};pvbo:{1}", Bus[gv.bus_no - gv.pv_no].g0, Bus[gv.bus_no - gv.pv_no].b0);
                            gv.Node_new.Add(gv.bus_no - gv.pv_no);
                            gv.Node_old[gv.bus_no - gv.pv_no] = i;
                            break;

                        case '3': //平衡节点
                            Bus[gv.bus_no].type = 3;
                            Bus[gv.bus_no].p = 0;
                            Bus[gv.bus_no].q = 0;
                            Bus[gv.bus_no].bus_serino = i;///add myself
                            Bus[gv.bus_no].vol = double.Parse(((string)tem_list[6]));
                            Bus[gv.bus_no].g0 = double.Parse((string)tem_list[16]);
                            Bus[gv.bus_no].b0 = double.Parse((string)tem_list[17]);
                            gv.Node_new.Add(gv.bus_no);
                            gv.Node_old[gv.bus_no] = i;
                            break;

                    }
                    line_a.Initialize();
                    tem_list.Clear();

                }

               
                /////////////读线路信息
                //读入第二部分：支路信息
                sr.ReadLine();
                line = sr.ReadLine();
                if (line != null)
                    line_a = line.Split(' ');
                for (int u = 0; u < line_a.Length; u++)
                {
                    if (line_a[u].Length != 0)
                    {
                        tem_list.Add(line_a[u]);
                    }
                }

                gv.branch_no = int.Parse((string)(tem_list[3]));
                line_a.Initialize();
                tem_list.Clear();
                Branch_type[] Branch = new Branch_type[gv.branch_no + 1];
                for (int i = 0; i <= gv.branch_no; i++)
                    Branch[i] = new Branch_type();

                string[] line_array = new string[gv.branch_no + 1];
                int maxlenght = 0;
                int indx = 1;

                ///////////////////////////////////////////////////////获取最大长度．
                for (int i = 1; i <= gv.branch_no; i++)
                {
                    line_array[i] = sr.ReadLine();
                    line = line_array[i];
                    line_a = line.Split(' ');
                    for (int u = 0; u < line_a.Length; u++)
                    {
                        if (line_a[u].Length != 0)
                        {
                            tem_list.Add(line_a[u]);
                        }
                    }
                    if (tem_list.Count > maxlenght && indx < 3)
                    {
                        maxlenght = tem_list.Count;
                        indx++;
                    }
                    tem_list.Clear();
                    line_a.Initialize();

                }
                for (int i = 1; i <= gv.branch_no; i++)
                {

                    line_a = line_array[i].Split(' ');
                    for (int u = 0; u < line_a.Length; u++)
                    {
                        if (line_a[u].Length != 0)
                        {
                            tem_list.Add(line_a[u]);
                        }
                    }

                    Branch[i].from = int.Parse((string)(tem_list[0]));
                    Branch[i].to = int.Parse((string)(tem_list[1]));

                    if (tem_list.Count < maxlenght)
                    {
                        Branch[i].r = double.Parse((string)(tem_list[4]));
                        Branch[i].x = double.Parse((string)(tem_list[5]));
                        Branch[i].b = double.Parse((string)(tem_list[6]));
                        Branch[i].k = double.Parse((string)(tem_list[12]));

                    }
                    else
                    {
                        Branch[i].r = double.Parse((string)(tem_list[5]));
                        Branch[i].x = double.Parse((string)(tem_list[6]));
                        Branch[i].b = double.Parse((string)(tem_list[7]));
                        Branch[i].k = double.Parse((string)(tem_list[13]));

                    }
                    Branch[i].start = (int)(gv.Node_new[Branch[i].from]);//节点编号重排序
                    Branch[i].end = (int)(gv.Node_new[Branch[i].to]);
                    line_a.Initialize();
                    tem_list.Clear();
                }
                ///////////////////////////////打印支路信息
                //for (int i = 1; i < Branch.Length; i++)
                //{
                //    Console.WriteLine("Branch" +"\t"+ Branch[i].from + "\t" + Branch[i].to + "\t" + Branch[i].r + "\t" + Branch[i].x + "\t" + Branch[i].b + "\t" + Branch[i].k + "\t" + Branch[i].start + "\t" + Branch[i].end);
                //}
                ////////////////////////////////
                /////////////////////////把所有母线和支路信息装入ｇｖ中的ｂｕｓ和ｂｒａｎｃｈ
                Bus_type bus0 = new Bus_type();
                gv.Bus.Add(bus0);//////为了使索引号从１开始
                Branch_type branch0 = new Branch_type();
                gv.Branch.Add(branch0);
                for (int i = 1; i < Bus.Length; i++)
                {
                    //  Console.WriteLine("serio_no:{0}", Bus[i].bus_serino);
                    gv.Bus.Add(Bus[i]);
                }
                for (int i = 1; i < Branch.Length; i++)
                {
                    gv.Branch.Add(Branch[i]);
                }

            }

            /////////////////////////打印母线信息
            StreamWriter sw1 = new StreamWriter("bus.txt");
            StreamWriter sww = new StreamWriter("scut300.txt");
            StreamWriter sww1 = new StreamWriter("fujian.txt");
            StreamWriter sww2 = new StreamWriter("fujian2.txt");
            string lll;
            ArrayList ll_list = new ArrayList();
            for (int j = 1; j < gv.Bus.Count; j++)
            {
                Bus_type bb = (Bus_type)(gv.Bus[j]);
                // Console.WriteLine(bb.bus_serino);
                // Console.WriteLine("第{6}个的type:{0} vol:{1} p_load:{2} q_load:{3} g0:{4} b0:{5}", bb.type, bb.vol, bb.p, bb.q, bb.g0, bb.b0, j);
                lll = "";
                lll = string.Format("第{6}个的type:{0} vol:{1} p_load:{2} q_load:{3} g0:{4} b0:{5}", bb.type, bb.vol, bb.p, bb.q, bb.g0, bb.b0, bb.bus_serino);
                sw1.WriteLine(lll);
                lll = "";
                double vol = bb.vol;
                if (bb.type == 2)
                    vol = -vol;
                lll = string.Format(" " + gv.Node_old[j] + " " + bb.p + " " + bb.q + " " + vol);
                sww.WriteLine(lll);
                lll = "";
                if (bb.g0 != 0)
                {
                    lll = string.Format(" " + gv.Node_old[j] + " " + 1 / bb.g0);
                    sww1.WriteLine(lll);
                    lll = "";
                    int start_index = 0;
                    lll = string.Format(" " + start_index + " " + gv.Node_old[j] + " " + 0 + " " + 0 + " " + bb.b0);
                    ll_list.Add(lll);//////把并联电容的母线信息存起来

                }
            }
            sw1.Close();
            sww1.Close();
            // Console.WriteLine("qq" + gv.Node_new[22]);
            ////////////////////////////////////////////////////////////////////////////////////
            /////////////////////////////////打印支路信息
            StreamWriter sw2 = new StreamWriter("branch.txt");
            //  Console.WriteLine("gv.Branch.Count"+gv.Branch.Count);
            for (int i = 1; i < gv.Branch.Count; i++)
            {
                Branch_type yy = (Branch_type)(gv.Branch[i]);
                lll = "";
                lll = string.Format("Branch" + "\t" + yy.from + "\t" + yy.to + "\t" + yy.r + "\t" + yy.x + "\t" + yy.b + "\t" + yy.k + "\t" + yy.start + "\t" + yy.end);
                sw2.WriteLine(lll);
                lll = "";
                double b = yy.b;
                double k = yy.k;
                int from = yy.from;
                if (yy.k != 0 && yy.k != 1)
                {
                    b = k;
                    from = -from;
                }
                lll = string.Format(" " + from + " " + yy.to + " " + yy.r + " " + yy.x + " " + b);
                sww.WriteLine(lll);
                if (from == 257 || from == 300 || yy.to == 257 || yy.to == 300)
                    sww2.WriteLine(lll);
                // Console.WriteLine("Branch" + "\t" + yy.from + "\t" + yy.to + "\t" + yy.r + "\t" + yy.x + "\t" + yy.b + "\t" + yy.k + "\t" + yy.start + "\t" + yy.end);
            }
            Console.WriteLine("ll_list.Count" + ll_list.Count);
            for (int i = 0; i < ll_list.Count; i++)
            {
                lll = "";
                lll = (string)ll_list[i];
                sww.WriteLine(lll);
            }
            sw2.Close();
            sww.Close();
            sww2.Close();
            //////////////////////////////
            ///////////////////////打印节点号对应表
            //for (int tt = 1; tt <= gv.bus_no; tt++)
            //{
            //    Console.WriteLine("节点号{0}，对应号{1}", (int)gv.Node_new[tt], tt);
            //}
            /////////////////////////

        }

        /*形成节点导纳矩阵*/
        public static void nodeAdmittanceMatrix()
        {
            int bus_no = gv.bus_no;
            int i = gv.i;

            for (i = 0; i < (gv.bus_no * gv.bus_no); i++)
                gv.NAM.Add(zero);//初始化并使ＮＡＭ有着bus_no * bus_no的维数；

            for (i = 0; i < gv.bus_no; i++) //母线对地导纳部分
            {

                Complex temp = new Complex(((Bus_type)(gv.Bus[i + 1])).g0, ((Bus_type)(gv.Bus[i + 1])).b0);
                gv.NAM[i * gv.bus_no + i] = temp;

            }
            //int genzong = 0;
            for (i = 1; i <= gv.branch_no; i++)//支路部分
            {

                Complex temp1 = new Complex(((Branch_type)gv.Branch[i]).r, ((Branch_type)gv.Branch[i]).x);
                Complex temp2 = new Complex();
                temp2 = 1 / temp1;//支路阻抗转换为导纳
                //if (i == 1)
                //    Console.WriteLine("temp2:"+temp2.Real+"+j"+temp2.Image);
                //pi型电路对地电纳
                Complex susceptance = new Complex(0, ((Branch_type)gv.Branch[i]).b);
                //以支路起始节点编号寻访一维的节点导纳矩阵
                //diag表对角线上元素　non_diag表非对角线元素．
                int diagi = (((Branch_type)gv.Branch[i]).start - 1) * bus_no + ((Branch_type)gv.Branch[i]).start - 1;
                int diagj = (((Branch_type)gv.Branch[i]).end - 1) * bus_no + ((Branch_type)gv.Branch[i]).end - 1;
                int non_diagi = (((Branch_type)gv.Branch[i]).start - 1) * bus_no + ((Branch_type)gv.Branch[i]).end - 1;
                int non_diagj = (((Branch_type)gv.Branch[i]).end - 1) * bus_no + ((Branch_type)gv.Branch[i]).start - 1;


                if (((Branch_type)gv.Branch[i]).k == 0)//普通支路
                {

                    gv.NAM[diagi] = (Complex)gv.NAM[diagi] + temp2 + susceptance * 0.5;//线路两端都有电纳，并联后去０．５
                    gv.NAM[diagj] = (Complex)gv.NAM[diagj] + temp2 + susceptance * 0.5;
                    //if (non_diagi == 1)
                    //{
                    //    Console.WriteLine(((Complex)gv.NAM[non_diagi]).Real + "+j" + ((Complex)gv.NAM[non_diagi]).Image);
                    //    Console.WriteLine("temp2:" +(-temp2).Real + "+j" +(-temp2).Image);
                    //}
                    gv.NAM[non_diagi] = (Complex)gv.NAM[non_diagi] + (-temp2);
                    gv.NAM[non_diagj] = (Complex)gv.NAM[non_diagj] + (-temp2);

                    //if (i == 1)
                    //    Console.WriteLine(((Complex)gv.NAM[non_diagi]).Real+"+j"+((Complex)gv.NAM[non_diagi]).Image);
                }
                else//变压器支路
                {
                    // genzong++;

                    double temp3 = 1 / ((Branch_type)gv.Branch[i]).k, temp4 = temp3 * temp3;
                    ///变压器绕组阻抗一率折算到始端且ｋ＝始端／末端
                    gv.NAM[diagi] = (Complex)gv.NAM[diagi] + temp2 * temp4 + susceptance * 0.5;

                    gv.NAM[diagj] = (Complex)gv.NAM[diagj] + temp2 + susceptance * 0.5;

                    gv.NAM[non_diagi] = (Complex)gv.NAM[non_diagi] + (-temp2) * temp3;

                    gv.NAM[non_diagj] = (Complex)gv.NAM[non_diagj] + (-temp2) * temp3;

                    //    if (genzong == 3)
                    //        Console.WriteLine("4-12trans:diagi: "+diagi+"k:"+temp3+"temp2.image:"+temp2.Image+"temp2.Real:"+temp2.Real+"z/k^2:"+(temp2*temp4).Real+"+j"+(temp2*temp4).Image);
                }
            }

            
        }
        ////////////////////////////////Matrix
        /*用高斯消元法求解线性化的以J阵为系数的方程组*/
        public static void Gauss(ref double[] jacobian, ref double[] x)
        {
            int j, k, num = gv.bus_no + gv.pq_no;//num就是Ｊ阵的维数即ｐｖ和ｐｑ节点的ｐ方程数＋ｐｑ节点的ｑ方程数；
            int i = gv.i;
            double temp;
            //消去过程,上三角化与首1化
            for (i = 1; i < num; i++)
            {
                temp = jacobian[(i - 1) * num + i - 1];//j阵对角元

                for (j = i; j <= num; j++)
                    jacobian[(i - 1) * num + j - 1] /= temp;//每行将对角元素单位化

                for (j = i; j < (num - 1); j++)
                {
                    temp = jacobian[j * num + i - 1];//i行以下的每行首个元素
                    for (k = i; k <= num; k++)
                        jacobian[j * num + k - 1] += -temp * jacobian[(i - 1) * num + k - 1];
                }
            }
            //回代过程
            for (i = (num - 1); i > 0; i--)
            {
                temp = 0;
                for (j = (num - 1); j > i; j--)
                {
                    temp += jacobian[(i - 1) * num + j - 1] * x[j - 1];
                }
                x[i - 1] = jacobian[(i - 1) * num + num - 1] - temp;
            }
        }
        ////////////////////////////////Gausss
        /*牛顿法迭代核心*/
        public static void Newton()
        {
            int bus_no = gv.bus_no;
            int pq_no = gv.pq_no;

            double[] delta_p = new double[bus_no];
            double[] p = new double[bus_no];
            double[] delta_q = new double[bus_no];
            double[] q = new double[bus_no];
            double[] delta_v = new double[pq_no + 1];
            double[] delta_o = new double[bus_no];
            double[] jacobian = new double[(bus_no + pq_no - 1) * (bus_no + pq_no)];//是增广的j阵
            double[] x = new double[bus_no + pq_no - 1];
            //////////////////////////////////////////////////////////////////////
            //Console.WriteLine("volvolvolvol");
            //for (int g = 1; g<= bus_no; g++)
            //{
            //    Console.WriteLine("No"+g+"\t\t"+((Bus_type)gv.Bus[g]).vol);
            //}
            //Console.WriteLine("volvolvolvol");
            /////////////////////////////////////////////////////////////////////
            bool flag = true; //iter = 0;
            int iter = gv.iter;
            int i = gv.i;
            while (flag)
            {
                int j, temp, temp1;
                double theta = 0;
                /////////////////////打印失配量
                //Console.WriteLine("///打印失配量");

                for (i = 1; i < bus_no; i++) //形成失配量
                {
                    double temp_p = 0, temp_q = 0;

                    for (j = 1; j <= bus_no; j++)
                    {

                        temp = (i - 1) * bus_no + j - 1;
                        theta = ((Bus_type)gv.Bus[i]).ang - ((Bus_type)gv.Bus[j]).ang;
                        temp_p = temp_p + ((Bus_type)gv.Bus[j]).vol * (Complex.GetReal(((Complex)gv.NAM[temp])) * Math.Cos(theta) + Complex.GetImage(((Complex)gv.NAM[temp])) * Math.Sin(theta));
                        temp_q = temp_q + ((Bus_type)gv.Bus[j]).vol * (Complex.GetReal(((Complex)gv.NAM[temp])) * Math.Sin(theta) - Complex.GetImage(((Complex)gv.NAM[temp])) * Math.Cos(theta));
                        //if (i == 1)
                        //{
                        //    Console.WriteLine("NO.j" + j + "\t" + theta + "\t\t" + temp_p + "\t" + temp_q + "\t" + ((Bus_type)gv.Bus[j]).vol);
                        //    Console.WriteLine("ppppppppp");
                        //}
                    }

                    p[i] = ((Bus_type)gv.Bus[i]).vol * temp_p;
                    delta_p[i] = ((Bus_type)gv.Bus[i]).p - p[i];

                    q[i] = ((Bus_type)gv.Bus[i]).vol * temp_q;
                    delta_q[i] = ((Bus_type)gv.Bus[i]).q - q[i];

                    //  Console.WriteLine(p[i] + "\t" + delta_p[i] + "\t" + q[i] + "\t" + delta_q[i]);
                }
                //Console.WriteLine("///打印失配量");
                //形成增广雅可比矩阵的各元素
                for (i = 1; i < bus_no; i++)//H阵的元素
                {
                    jacobian[(i - 1) * (bus_no + pq_no) + bus_no + pq_no - 1] = -delta_p[i];
                    for (j = 1; j < bus_no; j++)
                    {
                        temp = (i - 1) * bus_no + j - 1;
                        temp1 = (i - 1) * (bus_no + pq_no) + j - 1;
                        theta = ((Bus_type)gv.Bus[i]).ang - ((Bus_type)gv.Bus[j]).ang;
                        if (i == j)
                            jacobian[temp1] = ((Bus_type)gv.Bus[i]).vol * ((Bus_type)gv.Bus[i]).vol * Complex.GetImage(((Complex)gv.NAM[temp])) + q[i];
                        else
                            jacobian[temp1] = -((Bus_type)gv.Bus[i]).vol * ((Bus_type)gv.Bus[j]).vol * (Complex.GetReal(((Complex)gv.NAM[temp])) * Math.Sin(theta) - Complex.GetImage(((Complex)gv.NAM[temp])) * Math.Cos(theta));

                    }
                }

                for (i = 1; i < bus_no; i++)//N阵的元素
                    for (j = 1; j <= pq_no; j++)
                    {
                        temp = (i - 1) * bus_no + j - 1;
                        temp1 = (i - 1) * (bus_no + pq_no) + bus_no - 1 + j - 1;
                        theta = ((Bus_type)gv.Bus[i]).ang - ((Bus_type)gv.Bus[j]).ang;
                        if (i == j)
                            jacobian[temp1] = -((Bus_type)gv.Bus[i]).vol * ((Bus_type)gv.Bus[i]).vol * Complex.GetReal(((Complex)gv.NAM[temp])) - p[i];
                        else
                            jacobian[temp1] = -((Bus_type)gv.Bus[i]).vol * ((Bus_type)gv.Bus[j]).vol * (Complex.GetReal(((Complex)gv.NAM[temp])) * Math.Cos(theta) + Complex.GetImage(((Complex)gv.NAM[temp])) * Math.Sin(theta));
                    }

                for (i = 1; i <= pq_no; i++)//K阵的元素
                {
                    jacobian[(bus_no - 1 + i - 1) * (bus_no + pq_no) + bus_no + pq_no - 1] = -delta_q[i];
                    for (j = 1; j < bus_no; j++)
                    {
                        temp = (i - 1) * bus_no + j - 1;
                        temp1 = (bus_no - 1 + i - 1) * (bus_no + pq_no) + j - 1;
                        theta = ((Bus_type)gv.Bus[i]).ang - ((Bus_type)gv.Bus[j]).ang;
                        if (i == j)
                            jacobian[temp1] = ((Bus_type)gv.Bus[i]).vol * ((Bus_type)gv.Bus[i]).vol * Complex.GetReal(((Complex)gv.NAM[temp])) - p[i];
                        else
                            jacobian[temp1] = ((Bus_type)gv.Bus[i]).vol * ((Bus_type)gv.Bus[j]).vol * (Complex.GetReal(((Complex)gv.NAM[temp])) * Math.Cos(theta) + Complex.GetImage(((Complex)gv.NAM[temp])) * Math.Sin(theta));
                    }
                }

                for (i = 1; i <= pq_no; i++)//L阵元素
                    for (j = 1; j <= pq_no; j++)
                    {
                        temp = (i - 1) * bus_no + j - 1;
                        temp1 = (bus_no - 1 + i - 1) * (bus_no + pq_no) + bus_no - 1 + j - 1;
                        theta = ((Bus_type)gv.Bus[i]).ang - ((Bus_type)gv.Bus[j]).ang;
                        if (i == j)
                            jacobian[temp1] = ((Bus_type)gv.Bus[i]).vol * ((Bus_type)gv.Bus[i]).vol * Complex.GetImage(((Complex)gv.NAM[temp])) - q[i];
                        else
                            jacobian[temp1] = -((Bus_type)gv.Bus[i]).vol * ((Bus_type)gv.Bus[j]).vol * (Complex.GetReal(((Complex)gv.NAM[temp])) * Math.Sin(theta) - Complex.GetImage(((Complex)gv.NAM[temp])) * Math.Cos(theta));
                    }
                ////////////////
                //Console.WriteLine("jacojjjjjjjjjj");
                //for (i = 0; i < bus_no * bus_no; i++)
                //{
                //    if (jacobian[i] != 0)
                //        Console.WriteLine("第"+i+"个"+jacobian[i]);
                //}
                /////////////////
                FlowFunction.Gauss(ref jacobian, ref x);//高斯消去法求解线性化的以J阵为系数的方程组

                for (i = 1; i < bus_no; i++)//修正电压相角
                {
                    delta_o[i] = x[i - 1];
                    ((Bus_type)gv.Bus[i]).ang = ((Bus_type)gv.Bus[i]).ang + delta_o[i];
                }

                for (i = 1; i <= pq_no; i++)//修正电压幅值
                {
                    delta_v[i] = ((Bus_type)gv.Bus[i]).vol * x[bus_no + i - 2];
                    ((Bus_type)gv.Bus[i]).vol = ((Bus_type)gv.Bus[i]).vol + delta_v[i];
                }

                flag = false;//以下判定收敛条件
                for (i = 1; i < (bus_no - 1); i++)
                {
                    if (flag == false)
                        flag = (flag || (Math.Abs(delta_o[i]) > 0.00001));
                    else
                        break;
                }
                for (i = 1; i <= pq_no; i++)
                {
                    if (flag == false)
                        flag = (flag || (Math.Abs(delta_v[i]) > 0.00001));
                    else
                        break;
                }
                iter++; if (iter > gv.max_iter) break;/////////这里只设置了最多迭代３０次可以自己重新设置．
            }



            //将节点编号还原
            double[] temp_vol = new double[bus_no + 1];
            double[] temp_ang = new double[bus_no + 1];
            for (i = 1; i <= bus_no; i++)
            {
                temp_vol[i] = ((Bus_type)gv.Bus[((int)gv.Node_new[i])]).vol;
                temp_ang[i] = ((Bus_type)gv.Bus[((int)gv.Node_new[i])]).ang;
            }
            for (i = 1; i <= bus_no; i++)
            {
                ((Bus_type)gv.Bus[i]).vol = temp_vol[i];
                ((Bus_type)gv.Bus[i]).ang = temp_ang[i];
            }
            gv.iter = iter;
            // Console.WriteLine("bus_no="+bus_no+"iter="+iter);

        }
        //////////////////////newton function
        /////////////////////////////
        /*输出函数，计算支路功率及支路损耗*/
        public static void output()
        {
            int i = gv.i;
            StreamWriter sw = new StreamWriter("output.txt");
            string sq = "为了显示方便结果都只保留了三位小数进行了四舍五入";
            sw.WriteLine(sq);
            sq = string.Format("bus_no={0}\titer={1}", gv.bus_no, gv.iter);
            sw.WriteLine(sq);
            sw.WriteLine();
            sw.WriteLine();
            sq = string.Format("bus_no\t\tvoltage\t\tangle(degree)");
            sw.WriteLine(sq);

            for (i = 1; i <= gv.bus_no; i++)
            {
                sq = string.Format("{0}\t\t{1:F3}\t\t{2:F3}", i, ((Bus_type)gv.Bus[i]).vol, ((Bus_type)gv.Bus[i]).ang * 180 / Math.PI);
                sw.WriteLine(sq);
            }
            sw.WriteLine();
            sw.WriteLine();
            sq = string.Format("branch_no\tfrom\t\tto\t\tpower(i->j)\t\tpower(j->i)\t\tpower(loss)");
            sw.WriteLine(sq);


            for (i = 1; i <= gv.branch_no; i++)
            {
                //i是首端
                //j是末端
                // Console.WriteLine(((Branch_type)gv.Branch[i]).from);
                Complex_polar Vi = new Complex_polar(((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).from]).vol, ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).from]).ang);
                Complex_polar Vj = new Complex_polar(((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).to]).vol, ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).to]).ang);
                Complex vi = Complex_polar.change(Vi);//极坐标转换为直角坐标
                Complex vj = Complex_polar.change(Vj);
                Complex y0ij = new Complex(((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).start]).g0, ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).start]).b0 + ((Branch_type)gv.Branch[i]).b * 0.5);//母线上的并联导纳
                Complex y0ji = new Complex(((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).end]).g0, ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).end]).b0 + ((Branch_type)gv.Branch[i]).b * 0.5);
                Complex zij = new Complex(((Branch_type)gv.Branch[i]).r, ((Branch_type)gv.Branch[i]).x);
                Complex yij = new Complex(); ;
                yij = 1 / zij;//支路导纳
                if (((Branch_type)gv.Branch[i]).k == 0)//普通支路
                {
                    ((Branch_type)gv.Branch[i]).Sij = ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).from]).vol * ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).from]).vol * Complex.conj(y0ij) + vi * Complex.conj(vi - vj) * Complex.conj(yij);
                    ((Branch_type)gv.Branch[i]).Sji = ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).to]).vol * ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).to]).vol * Complex.conj(y0ji) + vj * Complex.conj(vj - vi) * Complex.conj(yij);
                }
                else//变压器支路
                {
                    ((Branch_type)gv.Branch[i]).Sij = ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).from]).vol * ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).from]).vol * Complex.conj(yij) * (1 - ((Branch_type)gv.Branch[i]).k) / (((Branch_type)gv.Branch[i]).k * ((Branch_type)gv.Branch[i]).k) + vi * Complex.conj(vi - vj) * Complex.conj(yij) / ((Branch_type)gv.Branch[i]).k;
                    ((Branch_type)gv.Branch[i]).Sji = ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).to]).vol * ((Bus_type)gv.Bus[((Branch_type)gv.Branch[i]).to]).vol * Complex.conj(yij) * (1 - 1 / ((Branch_type)gv.Branch[i]).k) + vj * Complex.conj(vj - vi) * Complex.conj(yij) / ((Branch_type)gv.Branch[i]).k;
                }
                ((Branch_type)gv.Branch[i]).delta_S = ((Branch_type)gv.Branch[i]).Sij + ((Branch_type)gv.Branch[i]).Sji;//支路损耗


            }

            for (i = 1; i <= gv.branch_no; i++)
            {
                if (Complex.GetImage(((Branch_type)gv.Branch[i]).Sij) >= 0 && Complex.GetImage(((Branch_type)gv.Branch[i]).Sji) >= 0)
                    sq = string.Format("{0}\t\t{1}\t\t{2}\t\t{3:F3}+j{4:F3}\t\t{5:F3}+j{6:F3}+\t\t{7:F3}+j{8:F3}", i, ((Branch_type)gv.Branch[i]).from, ((Branch_type)gv.Branch[i]).to, Complex.GetReal(((Branch_type)gv.Branch[i]).Sij), Complex.GetImage(((Branch_type)gv.Branch[i]).Sij), Complex.GetReal(((Branch_type)gv.Branch[i]).Sji), Complex.GetImage(((Branch_type)gv.Branch[i]).Sji), Complex.GetReal(((Branch_type)gv.Branch[i]).delta_S), Complex.GetImage(((Branch_type)gv.Branch[i]).delta_S));
                if (Complex.GetImage(((Branch_type)gv.Branch[i]).Sij) >= 0 && Complex.GetImage(((Branch_type)gv.Branch[i]).Sji) < 0)
                    sq = string.Format("{0}\t\t{1}\t\t{2}\t\t{3:F3}+j{4:F3}\t\t{5:F3}-j{6:F3}\t\t{7:F3}+j{8:F3}", i, ((Branch_type)gv.Branch[i]).from, ((Branch_type)gv.Branch[i]).to, Complex.GetReal(((Branch_type)gv.Branch[i]).Sij), Complex.GetImage(((Branch_type)gv.Branch[i]).Sij), Complex.GetReal(((Branch_type)gv.Branch[i]).Sji), -Complex.GetImage(((Branch_type)gv.Branch[i]).Sji), Complex.GetReal(((Branch_type)gv.Branch[i]).delta_S), Complex.GetImage(((Branch_type)gv.Branch[i]).delta_S));
                if (Complex.GetImage(((Branch_type)gv.Branch[i]).Sij) < 0 && Complex.GetImage(((Branch_type)gv.Branch[i]).Sji) >= 0)
                    sq = string.Format("{0}\t\t{1}\t\t{2}\t\t{3:F3}-j{4:F3}\t\t{5:F3}+j{6:F3}\t\t{7:F3}+j{8:F3}", i, ((Branch_type)gv.Branch[i]).from, ((Branch_type)gv.Branch[i]).to, Complex.GetReal(((Branch_type)gv.Branch[i]).Sij), -Complex.GetImage(((Branch_type)gv.Branch[i]).Sij), Complex.GetReal(((Branch_type)gv.Branch[i]).Sji), Complex.GetImage(((Branch_type)gv.Branch[i]).Sji), Complex.GetReal(((Branch_type)gv.Branch[i]).delta_S), Complex.GetImage(((Branch_type)gv.Branch[i]).delta_S));
                if (Complex.GetImage(((Branch_type)gv.Branch[i]).Sij) < 0 && Complex.GetImage(((Branch_type)gv.Branch[i]).Sji) < 0)
                    sq = string.Format("{0}\t\t{1}\t\t{2}\t\t{3:F3}-j{4:F3}\t\t{5:F3}-j{6:F3}\t\t{7:F3}+j{8:F3}", i, ((Branch_type)gv.Branch[i]).from, ((Branch_type)gv.Branch[i]).to, Complex.GetReal(((Branch_type)gv.Branch[i]).Sij), -Complex.GetImage(((Branch_type)gv.Branch[i]).Sij), Complex.GetReal(((Branch_type)gv.Branch[i]).Sji), -Complex.GetImage(((Branch_type)gv.Branch[i]).Sji), Complex.GetReal(((Branch_type)gv.Branch[i]).delta_S), Complex.GetImage(((Branch_type)gv.Branch[i]).delta_S));

                sw.WriteLine(sq);
            }

            sw.Close();

        }
    }

    class Program
    {
        static void Main(string[] args)
        {

            //int int_iter =0;
            // Console.WriteLine("请输入一个正整数作为最大迭代次数，默认为５０");
            // int_iter=Console.Read();
            //if (int_iter !=0)
            //   FlowFunction.gv.max_iter = int_iter;

            string FILE_NAME = "ieee300.txt";
            //FILE_NAME = Console.ReadLine();
            //FILE_NAME = FILE_NAME.Trim();
            DateTime dt1 = DateTime.Now;
            FlowFunction.Input(FILE_NAME);
            FlowFunction.nodeAdmittanceMatrix();
            FlowFunction.Newton();
            FlowFunction.output();
            DateTime dt2 = DateTime.Now;
            double dt = (dt2 - dt1).TotalSeconds;
            Console.WriteLine("iter=" + FlowFunction.gv.iter);
            Console.WriteLine("bus_no=" + FlowFunction.gv.bus_no);
            Console.WriteLine("在ｄｅｂｕｇ文件夹的ｏｕｔｐｕｔ．ｔｘｔ查看结果！");
            int jp;
            jp = Console.Read();
            Console.ReadLine();


        }
    }
}