C code for the space vector modulation algorithm

Question

I have been developing a C code for the space vector modulation of the three phase voltage source inverter. I have studied the SDK109 (page 93-109) and Digital Control in Power Electronics (Simone Buso, Paolo Mattavelli) and based on that I have attempted to write first version of the modulator. The code is intended to be used in a microcontroller but my idea was to first verify the behavior of my code in a simulation. So I have developed below given simulation in the Scilab/Xcos (free of charge equivalent of Matlab/Simulink).

The simulation consists of below given blocks:

1. Reference space vector generator

Generates reference space vector for the modulator. The magnitude of the reference space vector is equal to the radius of the inscribed circle to the hexagon normalized by the dc link voltage i.e. \$\frac{1}{\sqrt{3}}\$. The reference space vector rotates with angular velocity \$\omega=2\cdot\pi\cdot 60\,\mathrm{rad}\cdot\mathrm{s^{-1}}\$

2. SVM

Modulator itself i.e. block containing my C code of the space vector modulation algorithm. The C code inside

#include "scicos_block4.h"

#define U1  ((SCSREAL_COP *)GetRealInPortPtrs(block, 1))
#define U2  ((SCSREAL_COP *)GetRealInPortPtrs(block, 2))

#define Y1  ((SCSREAL_COP *)GetRealOutPortPtrs(block, 1))
#define Y2  ((SCSREAL_COP *)GetRealOutPortPtrs(block, 2))
#define Y3  ((SCSREAL_COP *)GetRealOutPortPtrs(block, 3))
#define Y4  ((SCSREAL_COP *)GetRealOutPortPtrs(block, 4))
#define Y5  ((SCSREAL_COP *)GetRealOutPortPtrs(block, 5))
#define Y6  ((SCSREAL_COP *)GetRealOutPortPtrs(block, 6))

#define X  ((SCSREAL_COP *)GetState(block))
#define dX ((SCSREAL_COP *)GetDerState(block))
#define Xk ((SCSREAL_COP *)GetDstate(block))
#define W  ((SCSREAL_COP *)GetWorkPtrs(block))

int getSector(double refVector[2][1]);
void calcSwtchStatesRelOnTimes(int sector, double refVector[2][1], double *delta1, double *delta2, double *delta3);
void calcDutyCycles(int sector, double *delta1, double *delta2, double *delta3, double *d1, double *d3, double *d5);

void multMatrices(double *matrix01, int noRows01, int noColumns01, double *matrix02, int noRows02, int noColumns02, double *result);
double scalarProductOfVectors(double *vector01, double *vector02, int length);
void getRow(double *matrix, int noRows, int noColumns, int row, double *dest);
void getColumn(double *matrix, int noRows, int noColumns, int column, double *dest);

void SVM(scicos_block *block,int flag)
{
  
  if (flag == 4) 
  {
    /* init */ 
   
  } 
  else if(flag == 1) 
  {
    /* output computation */ 
    double us_ref[][1] =
    {
        {0.0}, 
        {0.0}
    };

    // sector
    int sector;

    double delta1, delta2, delta3;
    double d1, d3, d5;


    // READING THE BLOCK INPUTS
    // -----------------------------------------------------------------------------------------------------------------

    us_ref[0][0] = U1[0];
    us_ref[1][0] = U2[0];


    // DETERMINATION OF THE SECTOR IN HEXAGON
    // -----------------------------------------------------------------------------------------------------------------
    sector = getSector(us_ref);

    // CALCULATION OF RELATIVE ON-TIME OF THE BASIC SPACE VECTORS
    // -----------------------------------------------------------------------------------------------------------------
    calcSwtchStatesRelOnTimes(sector, us_ref, &delta1, &delta2, &delta3);

    // CALCULATION OF DUTY CYCLES FOR INDIVIDUAL PHASES
    // -----------------------------------------------------------------------------------------------------------------
    calcDutyCycles(sector, &delta1, &delta2, &delta3, &d1, &d3, &d5);

    // WRITING TO THE BLOCK OUTPUTS
    // -----------------------------------------------------------------------------------------------------------------
    Y1[0] = d1;
    Y2[0] = d3;
    Y3[0] = d5;
    Y4[0] = sector;
    Y5[0] = delta1;
    Y6[0] = delta2;
  
  }
  else  if (flag == 5) 
  {
    /* ending */ 
  }

}

int getSector(double refVector[2][1])
{
    // transform matrix for coordinate system 1 
    double M1[][2] =
    {
            {1.0000, -0.5774},
            {0.0000,  1.1547}
    };

    // transform matrix for coordinate system 2 
    double M2[][2] =
    {
            { 1.0000, 0.5774},
            {-1.0000, 0.5774}
    };

    // transform matrix for coordinate system 3 
    double M3[][2] =
    {
            { 0.0000,  1.1547},
            {-1.0000, -0.5774}
    };
    
    // projections
    double temp[][1] = 
    {
            {0.0}, 
            {0.0}
    };

    // projections into basis space vectors
    double z1x, z1y, z2x, z2y, z3x, z3y;
    
    // sector
    int sector;
    
    // projections calculation

    // coordinate system 1
    multMatrices(&M1[0][0], 2, 2, &refVector[0][0], 2, 1, &temp[0][0]);
    z1x = temp[0][0]; 
    z1y = temp[1][0];

    // coordinate system 2
    multMatrices(&M2[0][0], 2, 2, &refVector[0][0], 2, 1, &temp[0][0]);
    z2x = temp[0][0]; 
    z2y = temp[1][0];

    // coordinate system 3
    multMatrices(&M3[0][0], 2, 2, &refVector[0][0], 2, 1, &temp[0][0]);
    z3x = temp[0][0]; 
    z3y = temp[1][0];

    // sector determination
    if(z1x*z1y >= 0)
    {
            if(z1x > 0)
            {
                    sector = 1;
            }
            else
            {
                    sector = 4;
            }
    }
    else
    {
            if(z2x*z2y >= 0)
            {
                    if(z2x > 0)
                    {
                            sector = 2;
                    }
                    else
                    {
                            sector = 5;
                    }
            }
            else
            {
                    if(z3x > 0)
                    {
                            sector = 3;
                    }
                    else
                    {
                            sector = 6;
                    }
            }
    }
    
    return sector;
}

void calcSwtchStatesRelOnTimes(int sector, double refVector[2][1], double *delta1, double *delta2, double *delta3)
{
  
    // transformation matrices for relative on-time calculations based on given reference space vector

    // sector 1
    double D1[][2] =
    {
            {0.0000,  1.0000},
            {0.8660, -0.5000}
    };

    // sector 2
    double D2[][2] =
    {
            { 0.8660,  0.5000},
            {-0.8660,  0.5000}
    };

    // sector 3
    double D3[][2] =
    {
            {-0.8660, -0.5000},
            { 0.0000,  1.0000}
    };

    // sector 4
    double D4[][2] =
    {
            {-0.8660,  0.5000},
            { 0.0000, -1.0000}
    };

    // sector 5
    double D5[][2] =
    {
            { 0.8660, -0.5000},
            {-0.8660, -0.5000}
    };

    // sector 6
    double D6[][2] =
    {
            {0.0000, -1.0000},
            {0.8660,  0.5000}
    };

    double temp[][1] = 
    {
            {0.0}, 
            {0.0}
    };
    
    double tmp1, tmp2, tmp3;
    
    // active switching states
    if(sector == 1)
    {
        multMatrices(&D1[0][0], 2, 2, &refVector[0][0], 2, 1, &temp[0][0]);
        tmp1 = temp[0][0];
        tmp2 = temp[1][0];
    }
    else
    {
        if(sector == 2)
        {
            multMatrices(&D2[0][0], 2, 2, &refVector[0][0], 2, 1, &temp[0][0]);
            tmp1 = temp[0][0];
            tmp2 = temp[1][0];
        }
        else
        {
            if(sector == 3)
            {
                multMatrices(&D3[0][0], 2, 2, &refVector[0][0], 2, 1, &temp[0][0]);
                tmp1 = temp[0][0];
                tmp2 = temp[1][0];
            }
            else
            {
                if(sector == 4)
                { 
                    multMatrices(&D4[0][0], 2, 2, &refVector[0][0], 2, 1, &temp[0][0]);
                    tmp1 = temp[0][0];
                    tmp2 = temp[1][0];
                }   
                else
                {
                    if(sector == 5)
                    {
                        multMatrices(&D5[0][0], 2, 2, &refVector[0][0], 2, 1, &temp[0][0]);
                        tmp1 = temp[0][0];
                        tmp2 = temp[1][0];
                    }
                    else
                    {
                        multMatrices(&D6[0][0], 2, 2, &refVector[0][0], 2, 1, &temp[0][0]);
                        tmp1 = temp[0][0];
                        tmp2 = temp[1][0];
                    }
                }
            }   
        }
    }
    
    // the rest of the pwm period is spent in zero switching states
    tmp3 = 1 - tmp1 - tmp2;   
    
    *delta1 = tmp1;
    *delta2 = tmp2;
    *delta3 = tmp3;
}

void calcDutyCycles(int sector, double *delta1, double *delta2, double *delta3, double *d1, double *d3, double *d5)
{
    double tmp1, tmp3, tmp5;
    
    if(sector == 1)
    {
        tmp1 = *delta3/2 + *delta1 + *delta2;
        tmp3 = *delta3/2 + *delta1;
        tmp5 = *delta3/2;
    }
    else
    {
        if(sector == 2)
        {
            tmp1 = *delta3/2 + *delta1;
            tmp3 = *delta3/2 + *delta1 + *delta2;
            tmp5 = *delta3/2;
        }
        else
        {
            if(sector == 3)
            {
                tmp1 = *delta3/2;
                tmp3 = *delta3/2 + *delta1 + *delta2;
                tmp5 = *delta3/2 + *delta1;
            }
            else
            {
                if(sector == 4)
                {
                    tmp1 = *delta3/2;
                    tmp3 = *delta3/2 + *delta1;
                    tmp5 = *delta3/2 + *delta1 + *delta2;
                }
                else
                {
                    if(sector == 5)
                    {
                        tmp1 = *delta3/2 + *delta1;
                        tmp3 = *delta3/2;
                        tmp5 = *delta3/2 + *delta1 + *delta2;
                    }
                    else
                    {
                        tmp1 = *delta3/2 + *delta1 + *delta2;
                        tmp3 = *delta3/2;
                        tmp5 = *delta3/2 + *delta1;
                    }
                }
            }
        }
    }
    
    *d1 = tmp1;
    *d3 = tmp3;
    *d5 = tmp5;
}

void multMatrices(double *matrix01, int noRows01, int noColumns01, double *matrix02, int noRows02, int noColumns02, double *result)
{
    int row, column;
    double *rowVector;
    double *columnVector;
    
    rowVector = (double *)malloc(sizeof(double)*noColumns01);
    columnVector = (double *)malloc(sizeof(double)*noRows02);
    
    for(row = 0; row < noRows01; row++)
    {
        for(column = 0; column < noColumns02; column++)
        {
            getRow(matrix01, noRows01, noColumns01, row, rowVector);
            getColumn(matrix02, noRows02, noColumns02, column, columnVector);
           
            *(result + row*noColumns02 + column) = scalarProductOfVectors(rowVector, columnVector, noColumns01);
        }
    }
}

double scalarProductOfVectors(double *vector01, double *vector02, int length)
{
    double product = 0.0;
    int elem;
    
    for(elem = 0; elem < length; elem++)
    {
        product += *(vector01 + elem)*(*(vector02 + elem));
    }
       
    return product;
}

void getRow(double *matrix, int noRows, int noColumns, int row, double *dest)
{
    int elem;
    for(elem = 0; elem < noColumns; elem++)
    {
        *(dest + elem) = *(matrix + row*noColumns + elem);
    }
}

void getColumn(double *matrix, int noRows, int noColumns, int column, double *dest)
{
    int elem;
    for(elem = 0; elem < noRows; elem++)
    {
        *(dest + elem) = *(matrix + column + elem*noColumns);
    }
}

3. PWM

Model of the pwm periphery of a microcontroller. Based on passed duty cycles for the three upper transistors in the inverter creates firing signals for those transistors.

4. VSI

Model of the voltage source inverter. Based on given firing signals for the upper three transistors in the inverter calculates phase and line voltages at the output of the inverter.

The waveforms of the duty cycles, firing signals, phase and line voltages look like that

1. Duty cycles

2. Firing signals for upper transistors along with the triangular carrier signal

3. Phase voltages

4. Line voltages

It seemed to me that the algorithm basically works. So I have attempted to filter out the first harmonic of the phase and line voltages with usage of the moving average filters with window width of 32 samples.

1. First harmonics of phase voltages

2. First harmonics of line voltages

According to the literature the amplitude of the first harmonic of the phase voltages should be \$\frac{U_{dc}}{\sqrt{3}}\$ and the amplitude of the first harmonic of the line voltages should be \$U_{dc}\$. I have set the \$U_{dc}=\sqrt{6}\cdot 230\,\mathrm{V}\$ (which corresponds to supplying the dc link from three phase grid 230V/50Hz via three phase diode bridge rectifier). So I have expected that the amplitude of the first harmonic of the phase voltage will be \$\sqrt{2}\cdot 230 \doteq 325\,\mathrm{V}\$ and I have found approximately \$120\,\mathrm{V}\$. For the line voltage I have expected \$U_{dc}=\sqrt{6}\cdot 230\doteq 563\,\mathrm{V}\$ and I have found approximately \$210\,\mathrm{V}\$.

It is obvious that the first harmonics of both voltages are smaller than they should be and that's my problem. I don't know what is the cause of this discrepancy. Does anybody see any mistake in my simulation - especially in the code for the SVM? Thanks in advance for any remarks.