Skip to content
Snippets Groups Projects
WindStressPRM.cs 18.6 KiB
Newer Older
  • Learn to ignore specific revisions
  • using System;
    using System.Collections.Generic;
    
    
    namespace WindStressPRM
    
        /// <summary>
        /// main calculations class
        /// </summary>
    
        public class StressPowerChecker
    
    Антон Кудряшов's avatar
    Антон Кудряшов committed
        {
    
            /// <summary>
            /// Input Data
            /// </summary>
    
            private Input Input;
    
            /// <summary>
            /// Main function for power graph algorithm
    
            /// Основая функция для вычисления распределения питания по графу ЛЭП-подстанции
    
            public Output CheckPower(Input input)
    
                this.Input = input;
    
                //Calculate which lines are broken
    
                List<WindStressPRM.Powerline> prmBrokenLines = brokenPowerLinesAfterCheck();
    
                foreach(PowerStation pwstation in this.Input.PowerStations) { pwstation.GetAttachedLines(this.Input.PowerLines); }
    
                PreparingPowerItems();
    
                foreach (PowerStation pwstation in input.PowerStations)
    
                    if (pwstation.IsSource)
    
                    {
                        CheckPowerPointsForStation(pwstation);
    
                Output output = new Output();
    
                output.DisabledStations = new List<PowerStation>();
    
                output.DisabledLines = new List<Powerline>();
    
                output.DisabledLines = prmBrokenLines;
    
                output.SpectificCoordinates = CheckSpecificLines();
    
                foreach (PowerStation powerStation in input.PowerStations)
    
                    //stations of type pole can be disabled if the line is broken 
    
                    if (!powerStation.IsON && !(powerStation.Stationtype == PowerStation.StationType.Pole))
    
                        output.DisabledStations.Add(powerStation);
    
                return output;
    
            /// <summary>
            /// recursive search function for power graph
            /// </summary>
            /// <param name="sourcepoint"> powered station to search next powered station from</param>
    
            private void CheckPowerPointsForStation(PowerStation sourcepoint)
    
                if (!sourcepoint.IsON)
    
                    throw new Exception("CheckPowerPointsForStation is called from disabled sourcepoint");
    
                foreach (Powerline line in sourcepoint.LineList)
    
                    // if attached powerline doesn't have power on, is not broken and has same or less voltage class than line currently powering station
                    // (currentVolt parameter) we turn this line's power ON
                    if (!line.IsBroken && !line.IsON) //&& line.Voltage <= sourcepoint.CurrentVolt
    
                        line.IsON = true;
                        foreach (PowerStation powerStation in Input.PowerStations)
    
                            // search to the powerstation that is connected to sourcepoint by the line
    
                            if (powerStation.Identifier != sourcepoint.Identifier && (powerStation.Identifier == line.PointFromID || powerStation.Identifier == line.PointToID))
    
                                // if the sourcepoint is not a pole  - i.e. we know 
                                // it can redistribute power to its  lines
                                // we turn connected point ON if any of the  connected lines are powered
                                if (!(sourcepoint.Stationtype == PowerStation.StationType.Pole) ) {
                                   // if (powerStation.CurrentVolt <= line.Voltage)
                                    {
                                        powerStation.SetCurrentVolt();
                                        powerStation.IsON = true;
                                        CheckPowerPointsForStation(powerStation);
                                    }
    
                                    // if line is a pole we have to check if it's actually able to
                                    // get electricity to other points i.e. no connected lines are broken
    
                                    bool powerLineCheck = false;
    
                                    foreach (Powerline powerline in powerStation.LineList)
    
                                        if (powerline.IsBroken) {
    
                                            powerLineCheck = true; 
                                        }
    
                                    if (!powerLineCheck)
    
                                      //  if (powerStation.CurrentVolt <= line.Voltage)
                                        {
                                            powerStation.SetCurrentVolt();
                                            powerStation.IsON = true;
                                            CheckPowerPointsForStation(powerStation);
                                        }
    
                    //else we have broken line or already switched ON powerpoint
                }
    
    
            /// <summary>
            /// preparing powerstations and lists of lines for them to get a graph-like structure
            /// </summary>
    
            private void PreparingPowerItems()
    
            {
                //First we make sure that all the sources are ON
                //and all non sources are OFF
    
                foreach (PowerStation powerStation in Input.PowerStations)
    
                    if (powerStation.IsSource == true) {
                        powerStation.IsON = true; 
    
                        powerStation.IsON = false; 
    
                    powerStation.TurnOffAttachedLines();
                    powerStation.InitCurrentVolt();
    
    
            /// <summary>
            /// function to chec if powerline is broken by wind
            /// </summary>
            /// <returns>list of broken powerlines</returns>
    
            private List<Powerline> brokenPowerLinesAfterCheck()
    
                List<Powerline> brokenLines = new List<Powerline>();
    
                // actually there are curves in powerLines
    
                foreach (Powerline powerCurve in Input.PowerLines)
    
                {
                    // get coordinates list
    
                    IList<Coordinate> points = powerCurve.Coordinates;
    
    
                    List<bool> checkList = new List<bool>();
                    // cycle throw all points in line
                    for (int i = 1; i < points.Count; i++)
                    {
                        double x1 = points[i - 1].X;
                        double y1 = points[i - 1].Y;
                        double x2 = points[i].X;
                        double y2 = points[i].Y;
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                        bool result = linearLineIsBroken(points[i - 1], points[i], powerCurve.Height, powerCurve.Voltage, powerCurve.Year);
    
                        checkList.Add(result);
                    }
                    foreach (bool chkpnt in checkList)
                    {
                        if (chkpnt == true)
                        {
    
                            powerCurve.IsBroken = true;
    
                    if (powerCurve.IsBroken)
    
                        brokenLines.Add(powerCurve);
    
            /// <summary>
            /// Divides the line between two vertices if needed
            /// depending on dist_threshold 
            /// and checks if line is broken by wind in any of midpoints
            /// </summary>
            /// <param name="coord1">first vertice coordinate</param>
            /// <param name="coord2">second vertice </param>
            /// <param name="heightLine">height</param>
            /// <param name="power">power for climatology switches</param>
            /// <returns>true if line is broken false otherwise</returns>
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            private bool linearLineIsBroken(Coordinate coord1, Coordinate coord2, double heightLine, int power, int year)
    
            {
                double distance = Math.Sqrt((coord2.X - coord1.X) * (coord2.X - coord1.X) + (coord2.Y - coord1.Y) * (coord2.Y - coord1.Y));
    
                double distpropD = distance / Input.DistThreshold;
    
                List<Coordinate> pointlist = new List<Coordinate>();
    
                Coordinate midpoint = new Coordinate(0, 0);
    
                int distpropI = Convert.ToInt32(distpropD);
                if (distpropI > 1)
                {
                    double constXdiff = (coord2.X - coord1.X) / distpropI;
                    double constYdiff = (coord2.Y - coord1.Y) / distpropI;
    
                    Coordinate subCoord1 = new Coordinate(coord1.X, coord1.Y);
                    Coordinate subCoord2 = new Coordinate(coord1.X + constXdiff, coord1.Y + constXdiff);
    
                    for (int j = 0; j < distpropI; j++)
                    {
    
                            // move to next center, at the next interval
                            subCoord1 = new Coordinate(subCoord2.X, subCoord2.Y);
                            subCoord2 = new Coordinate(subCoord1.X + constXdiff, subCoord1.Y + constYdiff);
    
                        midpoint = new Coordinate((subCoord1.X + subCoord2.X) / 2, (subCoord1.Y + subCoord2.Y) / 2);
                        pointlist.Add(midpoint);
    
                    midpoint = new Coordinate( (coord1.X + coord2.X) / 2, (coord1.Y + coord2.Y) / 2);
    
                Func<ClimateCell, double> climateClosure = delegate(ClimateCell cell) {
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                    if (year >= 1987)
                    {
                        return cell.Wind25;
                    }
    
                    else if (power > 5 && power <= 330)
    
                    {
                        return cell.Wind10;
    
                    else if (power <= 5)
                    {
                        return cell.Wind5;
    
                    else
                    {
                        return cell.Wind15;
                    }
                };
    
                foreach (Coordinate coords in pointlist)
    
                    double uwind = interpol(coords, Input.PrognosticCells, delegate(PrognosticCell cell) { return cell.VelocityX; });
                    double vwind = interpol(coords, Input.PrognosticCells, delegate(PrognosticCell cell) { return cell.VelocityY; });
                    double climwind = interpol(coords, Input.ClimateCells, climateClosure);
    
                    if (Double.IsNaN(uwind) || Double.IsNaN(vwind) || Double.IsNaN(climwind)) 
                    {
                        // interpolation fail. we can't say everything about here
                        // discussion: also we can save these points for detail analysis
                        res = false;
                        continue;
                    }
    
                    double umod = Math.Sqrt(uwind * uwind + vwind * vwind); ;
                    double angleline = Math.Atan2((coord2.Y - coord1.Y), (coord2.X - coord1.X));
                    double anglewind = Math.Atan2(vwind, uwind) - angleline;
                    double sinwind = Math.Sin(anglewind);
                    double C_height = 1.0;
    
                    { //wind is too low 
    
                    }
                    else
                    {   // calculate prognostic and threshold windstress
                        double p1 = -0.00078501;
                        double p2 = 0.13431;
                        double p3 = -2.11112;
                        double p4 = 19.548;
                        double qpr = p1 * umod * umod * umod + p2 * umod * umod + p3 * umod + p4;
                        double qcl = p1 * climwind * climwind * climwind + p2 * climwind * climwind + p3 * climwind + p4;
    
                        double Ppr = qpr * C_height *Input.GustsScalingFactor* sinwind * sinwind;
    
                        double Pcl = qcl * C_height * 1.0;
                        if (Ppr >= Pcl)
                        {
                            // here line is broken
    
                    checkbool.Add(res);
                }
                bool result = false;
                foreach (bool lastcheck in checkbool)
                {
    
                    if (lastcheck == true) 
                    { 
                        result = true;
                    }
    
            private List<Coordinate> CheckSpecificLines()
            {
                List<Coordinate> list = new List<Coordinate>();
                Coordinate firstCoordinate = Input.ClimateCells.Origin;
                Index lastIndex = new Index( Input.ClimateCells.RowsCount() - 1, Input.ClimateCells.ColumnCount() - 1);
                Coordinate lastCoordinate = Input.ClimateCells.CellToProjection(lastIndex);
                for (double x = firstCoordinate.X; x <= lastCoordinate.X; 
                    x += Math.Sign(lastCoordinate.X - firstCoordinate.X)*Input.Distance35kVCheck)
                {
                    for (double y = firstCoordinate.Y; y <= lastCoordinate.Y;
                    y += Math.Sign(lastCoordinate.Y - firstCoordinate.Y) * Input.Distance35kVCheck)
                    {
                        Coordinate current = new Coordinate(x, y);
                        double wind = interpol<ClimateCell>(current, Input.ClimateCells, delegate(ClimateCell cell) { return cell.Wind10; });
                        double u = interpol<PrognosticCell>(current, Input.PrognosticCells, delegate(PrognosticCell cell) { return cell.VelocityX; });
                        double v = interpol<PrognosticCell>(current, Input.PrognosticCells, delegate(PrognosticCell cell) { return cell.VelocityY; });
    
                        if (Double.IsNaN(u) || Double.IsNaN(v) || Double.IsNaN(wind))
                        {
                            // interpolation fail. we can't say everything about here
                            // discussion: also we can save these points for detail analysis
                            continue;
                        }
    
    
                        if (umod * Input.GustsScalingFactor > wind && umod >10) {
    
            private double interpol<T>(Coordinate coords, Matrix<T> matrix, Func<T, double> valueGetter)
    
            {
                // select directions for projections
                const bool normalX = true;// true - East, false West
                const bool normalY = false;// true - North, false South
    
                Index rc = matrix.ProjectionToCell(coords);
                Coordinate center = matrix.CellToProjection(rc);
    
                double xDiff = coords.X - center.X;
                double yDiff = coords.Y - center.Y;
                //calculate second index
                int row2, col2;
                if ((xDiff >= 0 && normalX) || (!normalX && xDiff < 0))
                {
    
                    row2 = rc.Row >= matrix.RowsCount() - 1 ? rc.Row - 1 : rc.Row + 1;
    
                }
                else
                {
                    row2 = rc.Row > 0 ? rc.Row - 1 : rc.Row + 1;
                }
                if ((yDiff >= 0 && normalY) || (!normalY && yDiff < 0))
                {
    
                    col2 = rc.Col >= matrix.ColumnCount() - 1 ? rc.Col - 1 : rc.Col + 1;
    
                }
                else
                {
                    col2 = rc.Col > 0 ? rc.Col - 1 : rc.Col + 1;
                }
                // indexes and values at bounds
    
                Index rcBotLeft = new Index(Math.Min(row2, rc.Row), Math.Min(col2, rc.Col));
                Index rcBotRight = new Index(Math.Max(row2, rc.Row), Math.Min(col2, rc.Col));
                Index rcTopLeft = new Index(Math.Min(row2, rc.Row), Math.Max(col2, rc.Col));
                Index rcTopRight = new Index(Math.Max(row2, rc.Row), Math.Max(col2, rc.Col));
    
                double valBotLeft = valueGetter(matrix.Cells[rcBotLeft.Row, rcBotLeft.Col]);
                double valBotRight = valueGetter(matrix.Cells[rcBotRight.Row, rcBotRight.Col]);
                double valTopLeft =  valueGetter(matrix.Cells[rcTopLeft.Row, rcTopLeft.Col]);
                double valTopRight = valueGetter(matrix.Cells[rcTopRight.Row, rcTopRight.Col]);
                Coordinate origin = matrix.CellToProjection(rcBotLeft);
    
                
                bool testBotLeft = Double.IsNaN(valBotLeft);
                bool testBotRight = Double.IsNaN(valBotRight);
                bool testTopLeft = Double.IsNaN(valTopLeft);
                bool testTopRight = Double.IsNaN(valTopRight);
    
                if (testBotLeft || testBotRight || testTopLeft || testTopRight) 
                {
                    // tests indicates that value at test-cell is missed. 
                    if (testTopRight && testTopLeft && testBotLeft && testBotRight) 
                    {
    
                        // Current value in such a bad place, you need to fill these place with raster values
                        return Double.NaN;
    
                    }
                    int count = 0;
                    if (testBotLeft)
                    {
                        valBotLeft = 0;
                        count++;
                    }
                    if (testBotRight)
                    {
                        valBotRight = 0;
                        count++;
                    }
                    if (testTopLeft)
                    {
                        valTopLeft = 0;
                        count++;
                    }
                    if (testTopRight)
                    {
                        valTopRight = 0;
                        count++;
                    }
                    //of course there is count not 0;
                    if (count == 0) {
                        throw new Exception("Interpolation Logic Error");
                    }
                    double average = (valTopLeft + valTopRight + valBotLeft + valBotRight)/count;
                    if (testBotLeft)
                    {
                        valBotLeft = average;
                    }
                    if (testBotRight)
                    {
                        valBotRight = average;
                    }
                    if (testTopLeft)
                    {
                        valTopLeft = average;
                    }
                    if (testTopRight)
                    {
                        valTopRight = average;
                    }
                }
    
                double hx = matrix.Size.Width;
                double hy = matrix.Size.Height;
    
                // coefficients
                double px = (coords.X - origin.X) / hx;
                double py = (coords.Y - origin.Y) / hy;
                // inverse directions
                px *= normalX ? 1 : -1;
                py *= normalY ? 1 : -1;
                // interpolation
                double top = (1 - px) * valTopLeft + px * valTopRight;
                double bot = (1 - px) * valBotLeft + px * valBotRight;
                double rval = (1 - py) * bot + py * top;
                return rval;
            }