Untitled

#include <bits/stdc++.h>
using namespace std;


int main() {

    int n, flag=0;
    printf("Number of Bus Bar: ");
    scanf("%d", &n);


    complex <double> Y[n+1][n+1];
    complex <double> tolarance (10e-10, 10e-10);

    printf("Consider zero(0) if there is no value\n\n");

    for(int i=1; i<=n; i++) {
        for(int j=1; j<=n; j++) {

            if(i>j) {

                Y[i][j].real() = Y[j][i].real();
                Y[i][j].imag() = Y[j][i].imag();
                continue;
            }

            printf("\nFor Y%d%d\n", i, j);
            printf("Conductive Part <space> Susceptive Part: ");

            scanf("%lf %lf", &Y[i][j].real(), &Y[i][j].imag());
        }
    }



    complex <double> slack;
    double slackMagnetude, slackAngle;

    printf("\n\nEnter Slack Bus Voltage (Magnetude <space> Angle): ");
    scanf("%lf %lf", &slackMagnetude, &slackAngle);

    slack.real(slackMagnetude*cos((slackAngle*3.1416)/180));
    slack.imag(slackMagnetude*sin((slackAngle*3.1416)/180));





    int pq, pv;
    printf("\n\nNumber of PQ Bus: ");
    scanf("%d", &pq);
    printf("Number of PV Bus: ");
    scanf("%d", &pv);
    printf("\n\n");

    if(1+pq+pv != n) {
        printf("Wrong Parameter inserted\n");
        return 0;
    }





    complex <double> PQ[n+1]    = { 0 };                  //for P & Q in complex
    complex <double> VDelP[n+1] = { 0 }, VDelPre[n+1];    //for V & Delta in polar
    complex <double> V[n+1]     = { 0 };                  //for V in Complex
    complex <double> Q[n+1]     = { 0 };                  //for Q-max & Q-min


    for(int i=2; i<=n; i++) {
        printf("\nFor Bus %d:\n", i);
        for(int j=1; j<=3; j++) {

            if(j==1) {
                printf("P%d (in PU) = ", i);
                scanf("%lf", &PQ[i].real());
            }

            else if(j==2 && i<=pq+1) {
                printf("Q%d (in PU) = ", i);
                scanf("%lf", &PQ[i].imag());
            }

            else if(j==3){
                if (i>pq+1) {
                    printf("|V%d| (in PU) = ", i);
                    scanf("%lf", &VDelP[i].real());

                    printf("Q (min): ");             //Taking Q-min & Q-max
                    scanf("%lf", &Q[i].imag());
                    printf("Q (max): ");
                    scanf("%lf", &Q[i].real());
                }
                else {
                    VDelP[i].real() = 1;            //initial value 1+0j
                    V[i].real() = 1;
                }


            }

        }
    }



    for(int itr=1; itr<=7; itr++) {

        printf("\n\nIteration: %d\n\n", itr);

        for(int i=2; i<=pq+1; i++) {                                     //for PQ Bus

            complex <double> temp;
            temp = (conj(PQ[i])/conj(V[i])) - (Y[i][1] * slack);

                for(int k=2; k<=n; k++) {

                    if(k==i) {continue;}


                    if(k<=pq+1) {
                        temp = temp - (Y[i][k]*V[k]);
                    }
                    else {
                        V[k].real() = VDelP[k].real() * cos((VDelP[k].imag()*3.1416)/180);
                        V[k].imag() = VDelP[k].real() * sin((VDelP[k].imag()*3.1416)/180);
                        temp = temp - Y[i][k]*V[k];
                    }
                }
                V[i] = temp/Y[i][i];

                printf("%Bus Voltage (V): %lf %lfj\n", V[i].real(), V[i].imag());
            }





            for(int i=pq+2; i<=n; i++) {                           //for PV bus

                complex <double> temp_sum;
                temp_sum = Y[i][1] * slack;

                for(int k=2; k<=n; k++) {                          //Q (Reactive Power)

                    if(k<=pq+1) {
                        temp_sum = temp_sum + Y[i][k]*V[k];
                    }

                    else {
                        V[k].real() = VDelP[k].real() * cos((VDelP[k].imag()*3.1416)/180);
                        V[k].imag() = VDelP[k].real() * sin((VDelP[k].imag()*3.1416)/180);
                        temp_sum = temp_sum + Y[i][k]*V[k];
                    }
                }
                PQ[i].imag() = (-1)* imag(conj(V[i]) * temp_sum);

                printf("Reactive Power (Q): %lf\n", PQ[i].imag());





                if(Q[i].imag() < PQ[i].imag() && PQ[i].imag() < Q[i].real()) {   // Delta (Angle)

                    complex <double> temp;
                    temp = (conj(PQ[i])/conj(V[i])) - (Y[i][1] * slack);

                    for(int k=2; k<=n; k++) {

                        if(k==i) {continue;}


                        if(k<=pq+1) {
                            temp = temp - (Y[i][k]*V[k]);
                        }
                    }
                    VDelP[i].imag() = (arg(temp/Y[i][i]));

                    printf("Voltage Angle (Delta): %lf\n\n", VDelP[i].imag()*180/3.1416);
                }

                else {                                                        //V (Magnitude + Angle)

                    if(PQ[i].imag() <= Q[i].imag()) {
                        PQ[i].imag() = Q[i].imag();
                    }
                    else if (PQ[i].imag() >= Q[i].real()) {
                        PQ[i].imag() = Q[i].imag();
                    }


                    complex <double> temp;
                    temp = (conj(PQ[i])/conj(V[i])) - (Y[i][1] * slack);

                    for(int k=2; k<=n; k++) {

                        if(k==i) {continue;}


                        if(k<=pq+1) {
                            temp = temp - (Y[i][k]*V[k]);
                        }
                        else {
                            V[k].real() = VDelP[k].real() * cos((VDelP[k].imag()*3.1416)/180);
                            V[k].imag() = VDelP[k].real() * sin((VDelP[k].imag()*3.1416)/180);
                            temp = temp - Y[i][k]*V[k];
                        }
                }
                V[i] = temp/Y[i][i];

                printf("%Bus Voltage (V): %lf %lfj\n", V[i].real(), V[i].imag());


                }
            }

            for(int t=0; t<n; t++) {
                if((abs(VDelP[t].real()-VDelPre[t].real()) <= tolarance.real())  &&  (abs(VDelP[t].imag()-VDelPre[t].imag()) <= tolarance.imag())) {
                    flag = 1;
                    continue;
                }
                flag = 0;
            }


            copy (VDelP, VDelP+n+1, VDelPre);

        }
}