How to calculate the integral and derivative term and integral and derivative time in a pid controller?

Question

Declaration:

I have asked same question here, I feel the community is big here, thus asked again. If there is a feature to link the questions, I would do that too:

https://engineering.stackexchange.com/questions/55184/how-to-calculate-the-integral-and-derivative-term-and-integral-and-derivative-ti

I am designing a PID controller.

1. While calculating integral controller, should you calculate mean of all the past errors or only calculate mean of past n samples? In the case of latter, how to pick n?

2. While calculating derivative controller, should you just fit a straight line to past m (for eg, past 5 points) points or take derivative of all the past errors and pass through a low pass filter (LPF)? In the case of latter, can anyone suggest a fast way to do that in C++ (derivative + LPF time not more than 0.3 ms).

3. About integral (Ti) and derivative (Td) time in the Standard form, How to choose the numbers? I am working on arresting atmospheric seeing related oscillations. What it means is when any astronomical object is imaged, due to atmosphere dynamics (random changes in refractive index of many stacked layers), the image position on the focal plane oscillates. These oscillations have typical frequencies of 500 Hz and in good conditions have smaller amplitudes about 1 arcsecond and bad conditions up to 4-6 arcsecond. The frame rate of my pid controller is 655 fps (655 Hz), which is faster than atmospheric seeing, but the question is how to decide the Ti and Td of pid controller in standard form?

Follows my current implementation of PID controller:

Here Kp, Ki, Kd are global variables.

Some other global variables are:

Ni - Number of points to compute integral error, if <0 means all past points, else past Ni points

Nd - Fits a straight line on past Nd points, the derivative is the slope of the fitted line.

XTemp and YTemp are image shifts of the current image with respect to the reference image, thus the correction will be negative of this.

limMinInt, limMaxInt are the min/max values integral controller can have, but neither I am not so sure of having this logic, nor I have any reasonable criteria to fix these values.

limMin, limMax are the min/max values of the output of the PID, however like above, do not know how to fix these values.

correctionShiftX, correctionShiftY are the values I want to find using the PID logic

*correctionX, *correctionY are the values given to actuator using the conversion matrix.

The function parameters are mostly pointers, thus are preserved from one iteration to next.

Please suggest if you can suggest improvement in this logic. Currently this PID logic, I am only able to arrest for about 3 minutes.

I believe a good pid tuning is required which also includes a good implementation.

inline int getCorrection(double XTemp, double YTemp, deque<double> &integralErrorX, deque<double> &integralErrorY,
                  deque<double> &derivativeErrorX, deque<double> &derivativeErrorY, uint64_t *integralCounter, uint64_t *derivativeCounter,
                  double *correctionX, double *correctionY, uint64_t curr_count, uint64_t nbImagesReceived
) {
    uint64_t i;
    double toCorrectXTemp, toCorrectYTemp;
    double xMean =  (double) (Nd + 1) / (double) 2;
    double xixbaryiybarX = 0, xixbaryiybarY = 0;
    double xixbarsq = 0;
    double integralErrorVoltageX = 0;
    double integralErrorVoltageY = 0;
    double derivativeErrorVoltageX = 0;
    double derivativeErrorVoltageY = 0;

    toCorrectXTemp = -1 * XTemp;
    toCorrectYTemp = -1 * YTemp;

    double correctionShiftX, correctionShiftY;

    if (Ni > 0) {
        if (integralErrorX.size() == Ni) {
            integralErrorX.pop_front();
            integralErrorY.pop_front();
        }
        integralErrorX.push_back(toCorrectXTemp);
        integralErrorY.push_back(toCorrectYTemp);

        if (integralErrorX.size() < Ni) {
            integralErrorVoltageX = 0;
            integralErrorVoltageY = 0;
        } else {
            integralErrorVoltageX = accumulate(integralErrorX.end() - Ni, integralErrorX.end(), 0.0) / Ni;
            integralErrorVoltageY = accumulate(integralErrorY.end() - Ni, integralErrorY.end(), 0.0) / Ni;
        }
    }
    else {
        if (integralErrorX.empty()) {
            integralErrorX.push_back(toCorrectXTemp);
            integralErrorY.push_back(toCorrectYTemp);
            integralErrorVoltageX = 0;
            integralErrorVoltageY = 0;
        }
        else {
            double sumX = integralErrorX.front();
            double sumY = integralErrorY.front();
            integralErrorX.pop_front();
            integralErrorY.pop_front();
            sumX += toCorrectXTemp;
            sumY += toCorrectYTemp;
            integralErrorX.push_back(sumX);
            integralErrorY.push_back(sumY);
            integralErrorVoltageX = sumX / (curr_count + 1);
            integralErrorVoltageY = sumY / (curr_count + 1);
        }
    }

    if (derivativeErrorX.size() == Nd) {
        derivativeErrorX.pop_front();
        derivativeErrorY.pop_front();
    }
    derivativeErrorX.push_back(toCorrectXTemp);
    derivativeErrorY.push_back(toCorrectYTemp);

    if (derivativeErrorX.size() < Nd) {
        derivativeErrorVoltageX = 0;
        derivativeErrorVoltageY = 0;
    }
    else {
        double meanDerivativeErrorVoltageX = 0;
        double meanDerivativeErrorVoltageY = 0;

        meanDerivativeErrorVoltageX = accumulate(derivativeErrorX.end()-Nd, derivativeErrorX.end(), 0.0) / Nd;
        meanDerivativeErrorVoltageY = accumulate(derivativeErrorY.end()-Nd, derivativeErrorY.end(), 0.0) / Nd;

        i = 0;
        for (auto it=derivativeErrorX.cbegin(); it!=derivativeErrorX.cend(); it++) {
            xixbaryiybarX += (*it - meanDerivativeErrorVoltageX) * (i + 1 - xMean);
            xixbarsq += ((i + 1 - xMean) * (i + 1 - xMean));
            i += 1;
        }
        i = 0;
        for(auto it=derivativeErrorY.cbegin(); it!=derivativeErrorY.cend(); it++){
            xixbaryiybarY += (*it - meanDerivativeErrorVoltageY) * (i + 1 - xMean);
            i += 1;
        }

        derivativeErrorVoltageX = xixbaryiybarX / xixbarsq;
        derivativeErrorVoltageY = xixbaryiybarY / xixbarsq;
    }

    double integralContributionX = Ki * integralErrorVoltageX;
    double integralContributionY = Ki * integralErrorVoltageY;

    if (integralContributionX <= limMinInt) {
        integralContributionX = limMinInt;
    }
    if (integralContributionX >= limMaxInt) {
        integralContributionX = limMaxInt;
    }
    if (integralContributionY <= limMinInt) {
        integralContributionY = limMinInt;
    }
    if (integralContributionY >= limMaxInt) {
        integralContributionY = limMaxInt;
    }
    correctionShiftX = Kp * toCorrectXTemp + integralContributionX + Kd * derivativeErrorVoltageX;
    correctionShiftY = Kp * toCorrectYTemp + integralContributionX + Kd * derivativeErrorVoltageY;

    if (correctionShiftX <= limMin) {
        correctionShiftX = limMin;
    }
    if (correctionShiftX >= limMax) {
        correctionShiftX = limMax;
    }
    if (correctionShiftY <= limMin) {
        correctionShiftY = limMin;
    }
    if (correctionShiftY >= limMax) {
        correctionShiftY = limMax;
    }

    *correctionX = A00*correctionShiftX + A01*correctionShiftY;
    *correctionY = A10*correctionShiftX + A11*correctionShiftY;

    sprintf(
        logString,
        "%ld, %ld, %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf, %lf",
        curr_count, nbImagesReceived, toCorrectXTemp, toCorrectYTemp,
        integralContributionX, integralContributionY, derivativeErrorVoltageX,
        derivativeErrorVoltageY, correctionShiftX, correctionShiftY,
        *correctionX, *correctionY
    );
    shift_uncorected_log(logString);
    return 0;
}