0
\$\begingroup\$

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;
}
\$\endgroup\$
  • \$\begingroup\$ I doubt there is a closed-form inverse function. Your best bet is probably to curve-fit a function to the numerical solution of the inverse function (easy to solve via binary search). It's relatively well behaved. Try something like 100 points over the 0-π range. Or simply interpolate a lookup table (or use a lookup table directly if your time resolution doesn't make that impractical). \$\endgroup\$ – Spehro Pefhany Dec 30 '18 at 3:27
3
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Very good. I wish you would re-write calc_alpha() to take a power level argument (0.0 to 1.0), and do away with Vrms and Vpk in this function. Like calc_alpha(float power). \$\endgroup\$ – mcu Dec 30 '18 at 12:20
  • \$\begingroup\$ Hmm, I am getting some wild output at certain input values. Can you try rms = 5.528 and pk = 10.0. \$\endgroup\$ – mcu Dec 30 '18 at 14:10
  • \$\begingroup\$ @mcu And there you go, the problem with Newton's method. There are bad starting points. With rms=5.528 and pk = 10, a starting point of alpha = 0.5 is bad. And we're dealing with cosine and sine which are extra bad because they are periodic. - If you change the starting point to 1, then it converges to "1.3943" which isn't even the right answer. Cool, who would've thought. - I thought that it would converge over a larger range. Oh well. \$\endgroup\$ – Harry Svensson Dec 30 '18 at 15:04
  • \$\begingroup\$ Maybe you can add a condition. If the equation diverges, instead of converging, try a different starting point. \$\endgroup\$ – mcu Dec 30 '18 at 15:08
  • \$\begingroup\$ @mcu I tried a couple and it always converged to 1.3943. - Not sure why it converges to an incorrect value. \$\endgroup\$ – Harry Svensson Dec 30 '18 at 15:10
2
\$\begingroup\$

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).

\$\endgroup\$
  • \$\begingroup\$ You are never using more than 5 digits, so you could use a float instead of a double for the same precision and also halving the space in memory. A float got 23 bits for the mantissa => \$log_{10}(2)×23\approx6.92\$. And 5 is less than 6.92. \$\endgroup\$ – Harry Svensson Dec 30 '18 at 10:01
  • \$\begingroup\$ A lookup table inside a microcontroller is indeed the end product of this task. But my question was more about how to arrive at this table in an elegant way. Can I ask you how you produced yours? \$\endgroup\$ – mcu Dec 30 '18 at 14:17
  • \$\begingroup\$ @HarrySvensson Indeed, in a real application the LUT would almost surely be converted to integers (eg. microseconds delay) \$\endgroup\$ – Spehro Pefhany Dec 30 '18 at 16:59
  • 1
    \$\begingroup\$ @SpehroPefhany Ahh, clever tricks. I underestimated your answer. \$\endgroup\$ – Harry Svensson Dec 30 '18 at 17:27
  • \$\begingroup\$ That will work. I have written something similar, only I was adjusting the angle by an ever decreasing amount to arrive at solution. Successive approximation. But I really wanted to have a real formula. Harry Svensson figured it out. Well, it is closer to being a formula I think. \$\endgroup\$ – mcu Dec 30 '18 at 17:32
1
\$\begingroup\$

I think this might work.

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

enter image description here

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

enter image description here enter image description here

for x = P

enter image description here

enter image description here

Well I can curve fit to your results.

\$\endgroup\$
  • \$\begingroup\$ Hi Tony. I added my comment to the main post. It is too long to post here. \$\endgroup\$ – mcu Dec 30 '18 at 5:29
  • \$\begingroup\$ Added 4th and 6th order polynomials instead.. What accuracy do you expect? Your results at 30 , 150 degrees indicate only 2.9% power which seems low. \$\endgroup\$ – Tony Stewart Sunnyskyguy EE75 Dec 30 '18 at 6:38
  • \$\begingroup\$ actually 2.9% avg is close enough on my sim with some triac drop when on for 30 deg. \$\endgroup\$ – Tony Stewart Sunnyskyguy EE75 Dec 30 '18 at 7:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.