Calculate Angle for Triac Phase Control

Question

Triac-based phase control of the mains power into resistive load.

According to this post,

https://www.physicsforums.com/threads/how-to-calculate-rms-voltage-from-triac-phase-angle.572668/#post-3737599

RMS voltage can be calculated as:

Vrms = Vpeak * sqr((2*pi - 2*angle + sin(2*angle)) / 4*pi)

I need the inverse of this. I need to solve this for angle.

I need to find the angle required to produce a required power level (specified as a fraction of 0.0 to 1.0 of full power).

While I can find the value of the angle by brute force, using the equation above, it would be better to have an equation that calculates the angle directly. My math is too rusty though.

Keywords: phase-fired controller, phase angle control.

--------------------------------------

C++ code, based on math done by Harry Svensson below.

#include <iostream>
#include <cmath>

#define pi 3.1415926535

using namespace std;

float f(float power, float alpha){
    return 2 * pi - 2 * alpha + sin(2 * alpha) - 4 * pi*pow(sqrtf(power) / sqrtf(2), 2);
    //calculate f
}

float df(float alpha){
    return 2 * cos(2 * alpha) - 2;
    //calculate the derivative of f
}

float calc_alpha(float power){
    float alpha = 0.5;
    //starting point = 0.5
    float d, oldD = 1000.0
    bool again = true;

    while (again) {
        for (int i = 0; i < 20; ++i){
            d = f(power, alpha) / df(alpha);
            if (abs(d) > abs(oldD)) {
                // the equation is diverging; the starting point was bad;
                alpha = (rand() % 3000) * 0.001 + 0.005;
                oldD = 1000.0;
                break;
            }
            oldD = d;
            alpha = alpha - d;
            if ((d<0.001) && (d>-0.001)){ return alpha; }
            if (i == 19) again = false;
        }
    }
    return alpha;
}

float calculatePower(float alpha){
    return (2 * pi - 2 * alpha + sin(2 * alpha)) / (4 * pi) * 2;
}

int main(){
    for (int i = 0; i < 180; i += 5){
        float alpha = i / 180.0 * pi;
        float power = calculatePower(alpha);
        float calculatedAlpha = calc_alpha(power);
        cout << alpha << "\t" << power << "\t" << calculatedAlpha << "\n";
    }
    cout << "Press any key to finish.";
    cin.get();
    return 0;
}

Answer 1

You want to find \$\alpha\$ in the following equation: $$ V_{rms} = V_{pk}\sqrt{\frac{2\pi - 2\alpha + sin(2\alpha)}{4\pi}}$$

In other words, you want to find out when: $$ V_{pk}\sqrt{\frac{2\pi - 2\alpha + sin(2\alpha)}{4\pi}}-V_{rms} = 0$$

Then you have an entire wikipage dedicated to Root-finding algorithms for solving these kinds of problems.

As the second equation is right now, I wouldn't go with Newton's method because I'm too lazy to find the derivative of the second equation with respect to \$\alpha\$. I would instead use the Secant method. But! I have a quick eye and can see that if I just rewrite the equation a little bit I can find the derivative of that very easily and then use Newton's equation. I do prefer to use Newton's equation when I can, because it converges pretty fast, compared to other methods.

---

So before we continue, allow me to go back to the first equation, divide both sides by \$V_{pk}\$ and then square both sides. This way we won't have to deal with square roots which is relatively time consuming. Let's multiply both sides by \$4\pi\$ while we're at it so we don't have to do any division. The function you have after doing what I just said is this one:

$$ 4\pi\bigg(\frac{V_{rms}}{V_{pk}}\bigg)^2 = 2\pi - 2\alpha + sin(2\alpha)$$

And again, let's find when it reaches 0 so we can apply the Root-finding algorithms.

$$ 2\pi - 2\alpha + sin(2\alpha) - 4\pi\bigg(\frac{V_{rms}}{V_{pk}}\bigg)^2 = 0$$

And let's derive the equation above with respect to \$\alpha\$ so we can use Newton's method:

$$ 2cos(2\alpha)- 2 = 0$$

Now we're ready to start rolling with Newton's method!

$$ x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)} \rightarrow \alpha_{n+1} = \alpha_n - \frac{2\pi - 2\alpha + sin(2\alpha) - 4\pi\bigg(\frac{V_{rms}}{V_{pk}}\bigg)^2}{2cos(2\alpha)- 2} $$

And we're done (ish), we just need to implement it in code. The equation above will converge quadratically to the correct \$\alpha\$.

---

Here is an example of Newton's method in action written in C++:

#include <iostream>
#include <cmath>

#define pi 3.1415926535

using namespace std;

float f(float Vrms,float Vpk, float alpha){
    return 2*pi-2*alpha+sin(2*alpha)-4*pi*pow(Vrms/Vpk,2);
    //calculate f
}

float df(float alpha){
    return 2*cos(2*alpha)-2;
    //calculate the derivative of f
}

float calc_alpha(float Vrms,float Vpk){

    float alpha=0.5;
    //starting point = 0.5

    for(int i=0;i<10;++i){
        //it usually converges within 10 steps

        float d = f(Vrms,Vpk,alpha)/df(alpha); 
        //Calculate how much we need to change so we can leave early
        //if we decide that we are close enough.

        cout << "\nalpha = " << alpha;

        alpha = alpha - d;
        //update our alpha 

        if((d<0.001) && (d>-0.001)){ return alpha;}
        //leave if we're close enough
        //otherwise we'll leave when the for loop ends
    }

    return alpha;
}

int main(){
    float rms = 3.79318159;
    float pk = 5.4;
    float a = calc_alpha(rms,pk);
    cout << "\n\nVrms = " << rms << "\nVpk = " << pk << "\nalpha = " << a << "\n";

    return 0;
}

Compiling the code above spits out the following information:

alpha = 0.5
alpha = 0.417462
alpha = 0.40068

Vrms = 3.79318
Vpk = 5.4
alpha = 0.400001

Which is correct, since I used Vpk = 5.4 and Alpha = 0.4 to calculate Vrms to 3.79318. And as you can see, it converged close enough in only 3 steps, on the 4th step alpha would've been equal to 0.4, for other values it might converge in less or more steps.

---

Edit:

After OP made some fiddling with the code above, he converged to this code:

#include <iostream>
#include <cmath>

#define pi 3.1415926535

using namespace std;

float f(float power, float alpha){
    return 2 * pi - 2 * alpha + sin(2 * alpha) - 4 * pi*pow(sqrtf(power) / sqrtf(2), 2);
    //calculate f
}

float df(float alpha){
    return 2 * cos(2 * alpha) - 2;
    //calculate the derivative of f
}

float calc_alpha(float power){
    float alpha = 0.5;
    //starting point = 0.5
    float d, oldD = 1000.0
    bool again = true;

    while (again) {
        for (int i = 0; i < 20; ++i){
            d = f(power, alpha) / df(alpha);
            if (abs(d) > abs(oldD)) {
                // the equation is diverging; the starting point was bad;
                alpha = (rand() % 3000) * 0.001 + 0.005;
                oldD = 1000.0;
                break;
            }
            oldD = d;
            alpha = alpha - d;
            if ((d<0.001) && (d>-0.001)){ return alpha; }
            if (i == 19) again = false;
        }
    }
    return alpha;
}

float calculatePower(float alpha){
    return (2 * pi - 2 * alpha + sin(2 * alpha)) / (4 * pi) * 2;
}

int main(){
    for (int i = 0; i < 180; i += 5){
        float alpha = i / 180.0 * pi;
        float power = calculatePower(alpha);
        float calculatedAlpha = calc_alpha(power);
        cout << alpha << "\t" << power << "\t" << calculatedAlpha << "\n";
    }
    cout << "Press any key to finish.";
    cin.get();
    return 0;
}

Which solves his problem. So this is a 2-man answer with my answer as a base, and OP "perfecting" it. I've gotten permission to copy & paste his solution into my answer so everything is in one place.

Answer 2

Using the super-dumb lookup table method, and using a dumb binary search to find values, I get:

double angle_lut[101] ={
0.00000, 0.02478, 0.03936, 0.05161, 0.06258, 0.07268,
0.08215, 0.09114, 0.09973, 0.10800, 0.11600, 0.12376,
0.13132, 0.13870, 0.14592, 0.15300, 0.15996, 0.16680,
0.17354, 0.18019, 0.18675, 0.19324, 0.19966, 0.20602,
0.21231, 0.21856, 0.22476, 0.23091, 0.23703, 0.24311,
0.24915, 0.25517, 0.26116, 0.26713, 0.27308, 0.27902,
0.28494, 0.29084, 0.29674, 0.30263, 0.30852, 0.31440,
0.32028, 0.32617, 0.33206, 0.33796, 0.34387, 0.34979,
0.35572, 0.36166, 0.36763, 0.37361, 0.37962, 0.38566,
0.39172, 0.39781, 0.40394, 0.41010, 0.41630, 0.42254,
0.42883, 0.43516, 0.44155, 0.44799, 0.45449, 0.46106,
0.46769, 0.47439, 0.48118, 0.48804, 0.49500, 0.50205,
0.50920, 0.51647, 0.52384, 0.53135, 0.53899, 0.54679,
0.55474, 0.56286, 0.57117, 0.57969, 0.58844, 0.59744,
0.60671, 0.61629, 0.62621, 0.63652, 0.64727, 0.65852,
0.67036, 0.68289, 0.69624, 0.71060, 0.72622, 0.74348,
0.76297, 0.78574, 0.81390, 0.85332, 1.00000};   

If power is an int from 0 to 100%, angle in radians is \$\pi\cdot\$(1.0-angle_lut[power])

You can interpolate the above table (linearly, or with a spline) if you want finer power levels, but usually 1% is good enough.

Edit: I encourage you to implement the reverse evaluation yourself, but the key step in my binary search implementation is this:

  angle_max = pi; angle_min = 0;          
  for (i = 0; i< 20; i++)
     {
     angle = (angle_max + angle_min)/2.0; 
     if (rms(angle) < power) 
         angle_max = angle; 
     else 
         angle_min = angle; 
     }

Here rms() is a function that returns the power fraction in terms of angle.

Since the interval is halved for each iteration of the loop, 20 iterations gives you \$2^{20}\$ or about 6 decimal digits in terms of angle.

This is about the dumbest search anyone would use (@HarrySvensson approach is much more intelligent and elegant) but it will reliably give you 6 digits of result in 20 iterations of the target function rms() vs. 500,000 average iterations if you stepped through it from beginning to target. Since the function is monotonic (as it needs to be over the interval of interest for this approach to work) there are other improvements possible in both cases, but modern computers are so fast that it would cost more in programmer time than machine time it would save (and the calculation is only done once, then it becomes constants in the code most likely).

Answer 3

I think this might work.

\$θ = {Arccos (1 - 2P )} \$ for per unit power (0 to 1)

Using Conduction angle from right to left to match direction as it appears on time scale

for x = P

Well I can curve fit to your results.