runeconomy-nx/runeconomy-xk/src/main/java/com/gyee/runeconomy/model/PowerFittingALG.java

package com.gyee.runeconomy.model;

import com.gyee.runeconomy.service.WindDirection.Point;

import java.math.BigDecimal;
import java.util.*;

/**
 * 功率曲线拟合算法
 * 最小二乘法
 */
public class PowerFittingALG {

    /**
     * 功率曲线拟合
     *
     * @param arrX        风速数组
     * @param arrY        功率数据
     * @param length      点个数
     * @param dimension   维度
     * @param scale       精度  0.1  0.01
     * @return
     */
    public static List<Point> buildLine(double[] arrX, double[] arrY, int length, int dimension, double scale) {

        List<Point> points = new ArrayList<>();

        if (arrX.length != arrY.length || arrX.length < 3) {
            return points;
        }

        double minValue = arrY[0];
        double maxValue = arrY[arrY.length - 1];

        double min = 0;
        double max = 0;

        double[] coefficient = MultiLine(arrX, arrY, length, dimension);

        for (double i = arrX[0]; i <= arrX[arrX.length - 1]; i += scale) {
            Point point = new Point();
            point.setX(i);


            for (int j = 0; j < coefficient.length; j++) {
                if (j == 0) {
                    point.setY(coefficient[j] * Math.pow(point.getX(), j));
                } else {
                    double temp = coefficient[j] * Math.pow(point.getX(), j);
                    point.setY(point.getY() + temp);
                }

            }
            if (point.getY() < minValue) {
                point.setY(minValue);

            }
            if (point.getY() > maxValue) {
                point.setY(maxValue);
            }

            if (point.getY() < min) {
                min = point.getY();
            }
            if (point.getY() > max) {
                max = point.getY();
            }

            points.add(point);
        }
        Builder(points, min, max);
        return points;
    }

    private static void Builder(List<Point> points, double min, double max) {
        boolean b = false;
        for (int i = 0; i < points.size(); i++) {
            if (b) {
                points.get(i).setY(max);
            } else {
                if (max == points.get(i).getY()) {
                    b = true;
                }
            }

        }

        for (int i = points.size() - 1; i > -1; i--) {
            if (!b) {
                points.get(i).setY(min);
            } else {
                if (min == points.get(i).getY()) {
                    b = false;
                }
            }

        }
    }


    ///<summary>
    ///用最小二乘法拟合二元多次曲线
    ///</summary>
    ///<param name="arrX">已知点的x坐标集合</param>
    ///<param name="arrY">已知点的y坐标集合</param>
    ///<param name="length">已知点的个数</param>
    ///<param name="dimension">方程的最高次数</param>
    //二元多次线性方程拟合曲线
    private static double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension) {
        int n = dimension + 1;                  //dimension次方程需要求 dimension+1个 系数
        double[][] Guass = new double[n][n + 1];      //高斯矩阵 例如：y=a0+a1*x+a2*x*x
        for (int i = 0; i < n; i++) {
            int j;
            for (j = 0; j < n; j++) {
                Guass[i][j] = SumArr(arrX, j + i, length);
            }
            Guass[i][j] = SumArr(arrX, i, arrY, 1, length);
        }
        return ComputGauss(Guass, n);
    }

    //求数组的元素的n次方的和
    private static double SumArr(double[] arr, int n, int length) {
        double s = 0;
        for (int i = 0; i < length; i++) {
            if (arr[i] != 0 || n != 0)
                s = s + Math.pow(arr[i], n);
            else
                s = s + 1;
        }
        return s;
    }

    private static double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length) {
        double s = 0;
        for (int i = 0; i < length; i++) {
            if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0))
                s = s + Math.pow(arr1[i], n1) * Math.pow(arr2[i], n2);
            else
                s = s + 1;
        }
        return s;

    }

    private static double[] ComputGauss(double[][] Guass, int n) {
        int i, j;
        int k, m;
        double temp;
        double max;
        double s;
        double[] x = new double[n];
        for (i = 0; i < n; i++) x[i] = 0.0;//初始化

        for (j = 0; j < n; j++) {
            max = 0;
            k = j;
            for (i = j; i < n; i++) {
                if (Math.abs(Guass[i][j]) > max) {
                    max = Guass[i][j];
                    k = i;
                }
            }


            if (k != j) {
                for (m = j; m < n + 1; m++) {
                    temp = Guass[j][m];
                    Guass[j][m] = Guass[k][m];
                    Guass[k][m] = temp;
                }
            }
            if (0 == max) {
                // "此线性方程为奇异线性方程" 
                return x;
            }

            for (i = j + 1; i < n; i++) {
                s = Guass[i][j];
                for (m = j; m < n + 1; m++) {
                    Guass[i][m] = Guass[i][m] - Guass[j][m] * s / (Guass[j][j]);
                }
            }

        }//结束for (j=0;j<n;j++)

        for (i = n - 1; i >= 0; i--) {
            s = 0;
            for (j = i + 1; j < n; j++) {
                s = s + Guass[i][j] * x[j];
            }
            x[i] = (Guass[i][n] - s) / Guass[i][i];
        }
        return x;
    }//返回值是函数的系数


    /**
     * 推力系数 CP 值
     * @param sweptarea  扫风面积
     * @param line
     * @return
     */
    public static LineCurveFitting buildCp(Double sweptarea, LineCurveFitting line){
        List<Point> cpValue = new ArrayList<>();

        double kqmd = 1.225; //空气密度
        double max = 0;
        double cpAvg = 0;
        double sum1 = 0;

        for (int i = 0; i < line.getYLines().size(); i++)
        {
            Point point = line.getYLines().get(i);
            double speed = point.getX();
            double power = point.getY();
            double k = 0.5 * kqmd * Math.pow(speed, 3) * sweptarea;
            double result = 0;
            if (k != 0)
                result = power / k * 1000;
            //s5.Points.AddXY(speed, result);
            Point cppoint = new Point();
            cppoint.setX(speed);
            cppoint.setY(result);
            cpValue.add(cppoint);

            if (max < result)
            {
                max = result;
            }

            if (result > 0 && speed >= 3 && speed <= 25)
            {
                cpAvg += result;
                sum1++;

            }

        }
        if(sum1>0)
            cpAvg /= sum1;

        line.setCpValue(cpValue);
        line.setCpAvg(new BigDecimal(cpAvg).setScale(2, BigDecimal.ROUND_CEILING).doubleValue());

        return line;
    }





    /**
     * 曲线偏差率
     * @param points1 风速功率数组
     * @param points2
     * @param maxp
     * @return
     */
    public static double curveDeviationRatio(List<Point> points1, List<Point> points2, double maxp) {
        double result = -0;
        double pc = 0;
        if (points1 != null && points1.size() > 0 && points2 != null && points2.size() > 0)
        {
            double count = 0;
            double sum = 0;
            for (Point point : points1){
                Optional<Point> p = points2.stream().filter(it -> it.getX() == point.getX()).findFirst();
                if (p.isPresent()){
                    sum += Math.pow(point.getY() - p.get().getY(), 2);
                    count ++;
                    pc += point.getY() - p.get().getY();
                }
            }
            sum = Math.sqrt(sum);
            count = Math.sqrt(count);
            maxp = maxp * count;
            if (maxp != 0)
                result = sum / maxp * 100;
            if (pc < 0)
                result = 0 - result;
        }
        return result;
    }

    /**
     * 曲线偏差率  正负偏差
     * @param points1 风速功率数组
     * @param points2 风速功率数组(保证功率)
     * @param maxp  区间内的最大保证功率
     * @param mins 最小风速
     * @param maxs 最大风速
     * @return
     */
    public static double curveDeviationRatio2(List<Point> points1, List<Point> points2, double maxp, double mins, double maxs) {
        double result = -0;
        double pc = 0;
        if (points1 != null && points1.size() > 0 && points2 != null && points2.size() > 0)
        {
            double count = 0;
            double sum = 0;
            for (Point point : points1){
                Optional<Point> p = points2.stream().filter(it -> it.getX() == point.getX()).findFirst();
                if (p.isPresent() && p.get().getX() >= mins && p.get().getX() < maxs){
                    sum += Math.pow(point.getY() - p.get().getY(), 2);
                    count ++;
                    pc += point.getY() - p.get().getY();
                }
            }
            sum = Math.sqrt(sum);
            count = Math.sqrt(count);
            maxp = maxp * count;
            if (maxp != 0)
                result = sum / maxp * 100;
            if (pc < 0)
                result = 0 - result;
        }
        return result;
    }

}