WindStressPRM.cs

using System;
using System.Collections.Generic;

namespace WindStressPRM
{
    /// <summary>
    /// main calculations class
    /// </summary>
    public class StressPowerChecker
    {
        /// <summary>
        /// Input Data
        /// </summary>
        private Input Input;

        /// <summary>
        /// Main function for power graph algorithm
        /// Основая функция для вычисления распределения питания по графу ЛЭП-подстанции
        /// </summary>
        public Output CheckPower(Input input)
        {
            this.Input = input;
            //Calculate which lines are broken
            List<WindStressPRM.Powerline> prmBrokenLines = brokenPowerLinesAfterCheck();
            //get the graph
            PreparingPowerItems();
            //start from source points
            foreach (PowerStation pwstation in input.PowerStations)
            {
                if (pwstation.IsSource)
                {
                    CheckPowerPointsForStation(pwstation);
                }
            }
            //fill output
            Output output = new Output();
            output.DisabledStations = new List<PowerStation>();
            output.DisabledLines = prmBrokenLines;
            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");
            }
            // if the point is not a pole  - i.e. we know 
            // it can redistribute power within connected lines
            // we turn it ON if any of the  connected lines are powered
            foreach (Powerline line in sourcepoint.LineList)
            {
                if (!line.IsBroken && !line.IsON)
                {
                    line.IsON = true;
                    foreach (PowerStation powerStation in Input.PowerStations)
                    {
                        if (powerStation.Identifier != sourcepoint.Identifier && (powerStation.Identifier == line.PointFromID || powerStation.Identifier == line.PointToID))
                        {
                            if (!(sourcepoint.Stationtype == PowerStation.StationType.Pole)) {
                                powerStation.IsON = true;
                                CheckPowerPointsForStation(powerStation);
                            }
                            else {
                                // 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)
                                {
                                    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; 
                }
                else {
                    powerStation.IsON = false; 
                } 
                // for each power station we create a list of powerlines it is attached to
                List<Powerline> lines = new List<Powerline>();

                foreach (Powerline line in Input.PowerLines)
                {
                    //we also switch OFF all lines
                    line.IsON = false;
                    if (line.PointFromID == powerStation.Identifier || line.PointToID == powerStation.Identifier) {
                        lines.Add(line);
                    }
                }
                powerStation.LineList = lines;
            }
        }
        /// <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;
                    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);
                }
            }

            return brokenLines;
        }

        /// <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>
        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++)
                {
                    if (j > 0)
                    {
                        // 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);
                }
            }
            else {
                midpoint = new Coordinate( (coord1.X + coord2.X) / 2, (coord1.Y + coord2.Y) / 2);
                pointlist.Add(midpoint);
            }
            Func<ClimateCell, double> climateClosure = delegate(ClimateCell cell) {
                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;
                }
            };
            List<bool> checkbool = new List<bool>();
            foreach (Coordinate coords in pointlist)
            {
                bool res = false;
                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;
                if (umod < 20)
                { //wind is too low 
                    res = false;
                }
                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 * sinwind * sinwind;
                    double Pcl = qcl * C_height * 1.0;
                    if (Ppr >= Pcl)
                    {
                        // here line is broken
                        res = true;
                    }
                }
                checkbool.Add(res);
            }
            bool result = false;
            foreach (bool lastcheck in checkbool)
            {
                if (lastcheck == true) 
                { 
                    result = true;
                }
            }
            return result;
        }

        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;
                }
            }
            // sizes for cell
            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;
        }
    }
}