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
foreach(PowerStation pwstation in this.Input.PowerStations) { pwstation.GetAttachedLines(this.Input.PowerLines); }
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 = 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);
}
}
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)
{
// 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;
}
else {
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;
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 < 10)
{ //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 *Input.GustsScalingFactor* 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 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);
double x = firstCoordinate.X;
double y = firstCoordinate.Y;
while (Math.Sign(lastCoordinate.X - firstCoordinate.X) * (lastCoordinate.X - x) >= 0)
{
while (Math.Sign(lastCoordinate.Y - firstCoordinate.Y) * (lastCoordinate.Y - y) >= 0)
{
Coordinate current = new Coordinate(x, y);
y += Math.Sign(lastCoordinate.Y - firstCoordinate.Y) * Input.Distance35kVCheck;
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; });
double umod = Math.Sqrt(u * u + v * v);
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) {
list.Add(current);
}
}
x += Math.Sign(lastCoordinate.X - firstCoordinate.X) * Input.Distance35kVCheck;
}
return list;
}
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;
}
}
}