Alex Nazarovsky

MATLAB, DSP, Julia, Power quality, Engineering, Metrology

Installing EPICS Core System and CSS Studio on Red Hat Linux

— Install EPICS BASE —

Become superuser

1
[user@computer ~]# su

Setup HTTP and HTTPS proxy (If you are not using one - just omit these lines)

1
2
[root@computer ~]# export http_proxy="http://user:password@server:port/"
[root@computer ~]# export https_proxy="http://user:password@server:port/"

Download EPICS base (3.15.4)

1
[root@computer ~]# wget https://www.aps.anl.gov/epics/download/base/base-3.15.4.tar.gz

Download CSS (4.1.1)

1
[root@computer ~]# wget https://ics-web.sns.ornl.gov/css/updates/apps/sns-css-4.1.1-linux.gtk.x86_64.zip

Make EPICS base directory

1
2
[root@computer ~]# mkdir /opt/epics
[root@computer ~]# cd /opt/epics

Extract files from downloaded EPICS archive

1
[root@computer epics]# tar -zxvf ~/base-3.15.4.tar.gz 

Make symbolic link /opt/epics/base

1
2
[root@computer epics]# ln -s base-3.15.4/ base
[root@computer epics]# cd base

Set EPICS CONFIG parameters (These parameters are for 64 bit version)

1
2
3
[root@computer base]# export EPICS_BASE=/opt/epics/base
[root@computer base]# export EPICS_EXTENSIONS=/opt/epics/extensions
[root@computer base]# export EPICS_ARCH="linux-x86_64"

Install readline library and run ‘make’

1
2
3
[root@computer base]# yum install readline-devel
[root@computer base]# make
[root@computer base]# export PATH=${PATH}:/opt/epics/base/bin/linux-x86_64

— Install CSS Studio —

Unzip files and change owner of CSS folder to ‘user’

1
2
3
[root@computer ~]# cd ~/
[root@computer ~]# unzip sns-css-4.1.1-linux.gtk.x86_64.zip -d /home/user/CSS/
[root@computer ~]# chown -R user /home/user/CSS

How to Model ADC Quantization Effect With Known Maximum Voltage and Number of Bits and Apply It to a Signal in Matlab

Model ADC quantization effect
1
2
3
   Umax=1; % Set maximum voltage, for which max ADC code is reached
   bits=16; % ADC resolution 
   U=double(uencode(U,bits,Umax,'signed'))/(2^(bits-1))*(Umax);

ADC quantization demo
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
t=0:1/500:1-1/500;
U=sin(2*pi*5.00*t);

Umax=1; % Set maximum voltage, for which max ADC code is reached
bits=3; % ADC resolution 
Uq=double(uencode(U,bits,Umax,'signed'))/(2^(bits-1))*(Umax);

figure;
subplot 211
plot(t,U), hold on, stem(t,Uq);
xlabel('t, s') , ylabel('Magnitude, V');
legend('U','U_{quantized}');
subplot 212
plot(t,U-Uq);
legend('Quantization Error');
xlabel('t, s') , ylabel('Magnitude, V');

Python Script to Download Gost From protect.gost.ru as a Pdf File

Unfortunately the Russian database of the national standards (GOSTs) does not enable to get files as .pdf , but only as separate images (which are even not directly available). So we can’t even download the texts, it’s a pity. I’ve made a simple Python script to address this problem, because it’s far better to read and explore a printed copy. Script uses some libraries, you can install them by running pip install img2pdf pypdf2 requests It was tested in Python 2.7.

Usage:

  1. Find the standard you want to download from site protect.gost.ru
  2. Get the link, in the example it is http://protect.gost.ru/v.aspx?control=8&baseC=-1&page=0&month=-1&year=-1&search=&RegNum=1&DocOnPageCount=15&id=126445
  3. Replace the link in the script
  4. Run the script python gostdl.py
  5. In the same folder you will get separate pages of document in .pdf and .jpg format as well as merged document in document.pdf
Python 2.7 script ‘gostdl.py’ to download Gost from protect.gost.ru as a pdf file
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import requests,re,img2pdf
from PyPDF2 import PdfFileMerger, PdfFileReader

s = requests.Session()
response = s.get("http://protect.gost.ru/v.aspx?control=8&baseC=-1&page=0&month=-1&year=-1&search=&RegNum=1&DocOnPageCount=15&id=126445")
links = re.findall(r"\" href=\".*?pageK=(.*?-.*?-.*?-.*?)\".*?>.*?</a>\s+.*?<a style", response.content)
n=1
pdf_merger = PdfFileMerger()

for link in links:
    url="http://protect.gost.ru/image.ashx?page="+link
    fname=str(n).zfill(2)+".jpg"
    print "Downloading "+ fname + " url=" +url
    r = s.get(url)
    with open(fname, "wb") as jpg:
       jpg.write(r.content)
       jpg.flush()

    pdf_bytes = img2pdf.convert([fname])
    pdf_name = str(n).zfill(2)+".pdf"
    with open(pdf_name,"wb") as pdf:
       pdf.write(pdf_bytes)
    pdf_merger.append(PdfFileReader(file(pdf_name, 'rb')))
    n=n+1

pdf_merger.write("document.pdf")

Fixing Windows From BSOD (STOP 7B 0x0000007B) After a Motherboard Replacement

After the motherboard replacement I’ve got a blue screen of death (BSOD). This can happen after hard disk replacement also. I’ve used the recipe from here to recover from this failure.

  1. Download the script file ( mirror )
  2. Put the script file into the root of your USB Flash
  3. Boot to the Recovery Console (F8 > Repair Your Computer) or from a system repair/Windows 7 disc.
  4. Launch the command prompt.
  5. To determine the name of the USB Flash disk you can use wmic logicaldisk get name,description
  6. Change directory to the root of USB disk, e.g cd E:\
  7. Run the script cscript fix_7hdc.vbs /enable /search

How to Effectively Convert Books From DJVU to PDF

DJVU is a very good, although not very portable format for books and documents. PDF is less efficient in compression, but is far more popular and portable. You can easily find a free PDF reader for any system, for example for iPhone or iPad, but opening and reading DJVU is harder.

But you can always quiclky export DJVU to PDF and then use your favourite PDF reader program. On the iPad I use GoodReader and FoxitReader.

Here are the instructions:

1 - Install DJVULibre. It is available for Windows, Linux and MacOS.

2 - Use File -> Export As…

3 - Choose export settings

4 - Choose compression settings

5 - You can also convert books from the command line (good for batch procesing)

command line DJVU to PDF conversion (B & W mode good for most of the books)
1
ddjvu.exe -format=pdf -mode=black -quality=75 infile.djvu outfile.pdf 

Calculate Powers in Non-sinusoidal Conditions According to IEEE 1459

Definitions of electric power qualities under nonsinusoidal conditions are quite complex and hard to grasp. Annex B of the IEEE Std 1459-2010 shows how to calculate apparent power components in non-sinusoidal conditions. I have written a MATLAB script to illustrate this example. Unfortunately, there were typo errors in the standard, that’s why some of the values are different in my example.

Figure 1 — Instantaneous voltage and current

A load is supplied with a nonsinusoidal voltage

\[ v=v_1 + v_3 + v_5 +v_7 = \sqrt{2} \sum_{h=1,3,5,7}V_h \sin(h \omega t-\alpha_h) \]

has a nonsinusoidal current

\[ i=i_1 + i_3 + i_5 +i_7 = \sqrt{2} \sum_{h=1,3,5,7}I_h \sin(h \omega t-\beta_h) \]

(To simplify the explanations, the eventual presence of dc components was ignored.) In this case, the instantaneous power has 16 terms that can be separated in two groups

\[ p=vi=p_{hh}+p_{mn} \]

where

\[ p_{hh}=v_1 i_1 + v_3 i_3 +v_5 i_5 +v_7 i_7 = \sum_{h=1,3,5,7} v_h i_h \]

is the instantaneous power that contains only direct products (i.e., each component is the result of interaction of voltage and current harmonics of the same order).

\[ p_{mn}=v_1 (i_3 +i_5 +i_7) + v_3 (i_1 +i_5 +i_7) +v_5 (i_1 +i_3 +i_7) +v_7 (i_1+i_3+i_7) = \] \[ = \sum_{m=1,3,5,7} v_m \sum_{n =1,3,5,7 \atop n \neq m} i_n \]

is the instantaneous power that contains only cross products (i.e., each component is the result of interaction of voltage and current harmonics of different orders).

The direct products yield

\[ v_h i_h=\sqrt{2} \, V_h \sin (h \omega t - \alpha_h) \sqrt{2} \, I_h \sin (h \omega t - \beta_h) = \] \[ = P_h [1 - \cos(2 h \omega t - 2 \alpha_h)] - Q_h \sin (2 h \omega t - 2 \alpha_h) \]

where

\[ P_h = V_h I_h \cos (\theta_h)\] \[ Q_h =V_h I_h \sin (\theta_h) \]

are the harmonic active and reactive powers of order $ h $, respectively, and $ \theta_h = \beta_h - \alpha_h $ is the phase angle between the phasors $ V_h $ and $ I_h $.

The total active power is

\[ P = \sum_{h=1,3,5,7} P_h =P_1 +P_H \]

where

\[ P_1 = V_1 I_1 \cos(\theta_1) \]

is the fundamental (power-frequency) active power, and

\[ P_H=P_3+P_5+P_7 = \sum_{h \neq 1} P_h \]

is the total harmonic active power.

For each harmonic order, there is an apparent power of order h

\[ S_H=\sqrt{P_h^2 +Q_h^2} \]

The cross-products of the instantaneous powers have expressions as follows:

\[ v_m i_n=\sqrt{2} \, V_m \sin (m \omega t - \alpha_m) \sqrt{2} \, I_n \sin (n \omega t - \beta_n) = \] \[ = D_{mn} \left( \cos[(m-n) \omega t - \alpha_m + \beta_n ] + \cos[(m+n) \omega t - \alpha_m - \beta_n] \right) \]

where

\[ D_{mn}=V_m I_n \]

The total apparent power squared

\[ S^2 = V^2 I^2 =(V_1^2 +V_3^2 +V_5^2 +V_7^2)(I_1^2 +I_3^2 +I_5^2 +I_7^2) \]

may be separated in the same manner as the instantaneous power, in direct and the cross-products:

\[ S^2 = V_1^2 I_1^2 +V_3^2 I_3^2 + V_5^2 I_5^2 + V_7^2 I_7^2 + \] \[ + V_1^2 (I_3^2 +I_5^2 +I_7^2) + I_1^2 (V_3^2 +V_5^2 +V_7^2)+ \] \[ + V_3^2 I_5^2 + V_3^2 I_7^2 + V_5^2 I_3^2 + V_5^2 I_7^2 + V_7^2 I_3^2 + V_7^2 I_5^2 \]

or

\[ S^2 = S_1^2 +S_3^2 +S_5^2 +S_7^2 + \] \[ + D_I^2 +D_V^2 + \] \[ + D_{35}^2+ D_{37}^2 + D_{53}^2 + D_{57}^2 + D_{73}^2 + D_{75}^2 = S_1^2 + S_N^2 \]

where

\[ S_1^2 = P_1^2 + Q_1^2 \]

with $ S_1 $ , $ P_1 $ and $ Q_1 $ are the apparent, active, and reactive fundamental powers, and

\[ S_N^2 = D_I^2 + D_V^2 +S_H^2 \]

where

\[ D_I^2 = V_1^2 (I_3^2 +I_5^2 +I_7^2) \]

is the current distortion power,

\[ D_V^2 = I_1^2 (V_3^2 +V_5^2 +V_7^2) \]

is the voltage distortion power, and

\[ S_H^2 = S_3^2 +S_5^2 +S_7^2 + \] \[ + D_{35}^2+ D_{37}^2 + D_{53}^2 + D_{57}^2 + D_{73}^2 + D_{75}^2 = \] \[ = P_3^2 + P_5^2 + P_7^2 + Q_3^2 + Q_5^2 + Q_7^2 + \] \[ + D_{35}^2+ D_{37}^2 + D_{53}^2 + D_{57}^2 + D_{73}^2 + D_{75}^2 \qquad (1) \]

is the harmonic apparent power

If the load is supplied by a line with a resistance $ r $ the power loss in the line is

\[ \Delta P = r I^2 = \frac{r}{V^2} S^2 = \frac{r}{V^2} (S_1^2 + S_N^2) = \] \[ = \frac{r}{V^2} (P_1^2 +Q_1^2 +D_I^2 +D_V^2 +S_H^2) \qquad (2) \]

It is learned from this expression that every component of $ S $ contributes to the total power loss in the supplying system. This means that not only fundamental active and reactive powers cause losses but also the current and voltage distortion powers as well as the harmonic apparent power cause losses.

The following numerical example is meant to facilitate the understanding of the previous explanations:

The instantaneous voltages and currents are

\[ \begin{array}{c c} v_1=\sqrt{2} \, 100 \sin (\omega t - 0^{\circ} ) & i_1=\sqrt{2} \, 100 \sin (\omega t - 30^{\circ} ) \\ v_3=\sqrt{2} \, 8 \sin (\omega t - 70^{\circ} ) & i_3=\sqrt{2} \, 20 \sin (3 \omega t - 165^{\circ} ) \\ v_5=\sqrt{2} \, 5 \sin (\omega t + 140^{\circ} ) & i_5=\sqrt{2} \, 15 \sin (5 \omega t + 233^{\circ} ) \\ v_7=\sqrt{2} \, 7 \sin (\omega t + 20^{\circ} ) & i_7 =\sqrt{2} \, 10 \sin (7 \omega t + 288^{\circ} ) \\ \end{array} \]

The calculated active powers are summarized in Table 1.

Table 1 — Active powers

\[ \begin{array}{|c|c|c|c|c|c|} \hline P_1 (W) & P_3 (W) & P_5 (W) & P_7 (W) & P (W) & P_H (W) \\ \hline 8660.25 & -13.94 & -11.78 & -1.74 & 8632.79 & -27.47 \\ \hline \end{array} \]

The total harmonic active power $ P_H =P –27.47 W < 0 $ is supplied by the load and injected into the power system. This condition is typical for dominant nonlinear loads. The bulk of the active power is supplied to the load by the fundamental component $ P_1 = 8660.25 W$ .

The four reactive powers are given in Table 2.

Table 2 — Reactive powers

\[ \begin{array}{|c|c|c|c|} \hline Q_1 (var) & Q_3 (var) & Q_5 (var) & Q_7 (var) \\ \hline 5000.00 & 159.39 & -224.69 & 49.97 \\ \hline \end{array} \]

Of interest is the fact that $ Q_5 < 0$ , whereas other reactive powers are positive. If one incorrectly defines a total reactive power as the sum of the four reactive powers (in accordance with C. Budeanu’s definition):

\[ Q_B = Q_1 +Q_3 +Q_5 +Q_7 = 4984.67 \, var \]

and assumes that the supplying line has a resistance $ r= 1.0 \, \Omega $ and the load is supplied with an rms voltage $ V=240 V $, the power loss due to $ Q_B $ in line is

\[ \Delta P_B = \frac{r}{V^2}Q_B^2= \frac{1}{240^2}4984.67^2 = 431.37 \, W \]

According to the previous analysis, [see Equation (1) and Equation (2)], the correct way to find the corresponding power loss due to $ Q_1 $ , $ Q_3 $ , $ Q_4 $ , $ Q_7 $ is

\[ \Delta P = \frac{r}{V^2}(Q_1^2 +Q_3^2 +Q_5^2 +Q_7^2)= 435.39 \, W > P_B \]

The reactive power $ Q_5 $ , despite its negative value, contributes to the line losses in the same way as the positive reactive powers. The fact that harmonic reactive powers of different orders oscillate with different frequencies reinforces the conclusion that the reactive powers should not be added arithmetically (as recommended by Budeanu).

The cross-products that produce the distortion powers $ D_I $ and $ D_V $ are given in Table 3.

Table 3 — Distortion powers and their components

\[ \begin{array}{|c|c|c|c|} \hline D_{13} (var) & D_{15} (var) & D_{17} (var) & D_{I} (var) \\ \hline 2000.00 & 1500.00 & 1000.00 & 2692.58 \\ \hline D_{31} (var) & D_{51} (var) & D_{71} (var) & D_{V} (var) \\ \hline 800.00 & 1500.00 & 500.00 & 1772.00 \\ \hline \end{array} \]

Finally the remaining cross-products that belong to the harmonic apparent power are presented in Table 4.

Table 4 — Distortion harmonic powers

\[ \begin{array}{|c|c|c|c|c|c|} \hline D_{35} (var) & D_{37} (var) & D_{53} (var) & D_{57} (var) & D_{73} (var) & D_{75} (var) \\ \hline 120.00 & 80.00 & 300.00 & 150.00 & 100.00 & 75.00 \\ \hline \end{array} \]

The studied system has the rms voltage and current: $ V= 101.56 \, V $ and $ I= 103.56 \, A $ with the total harmonic distortions $ THD_V = 0.177 $ and $ THD_I = 0.269 $

The apparent power and its components are represented in the following tree:

Figure 2 — Apparent power and its components tree

The fundamental power factor (displacement power factor) is $ PF_1 = P_1 / S_1 =0.866 $, and the power factor is $ PF= P/S=0.821 $ . The dominant power components are $ P_1 $ and $ Q_1 $ . Due to relatively large distortion, $ S_N $ is found to be a significant portion of $ S $, and the current distortion power $ D_I $ is the dominant component of $ S_N $ .

Here is the MATLAB script that illustrates the previous material.

Example illustration to IEEE 1459-2010 Annex B
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
clear all; clc;
f_d=12800; % sampling frequency
F=50; % Mains frequency
t=0:1/f_d:1-1/f_d; % time
w=2*pi*F; % circular frequency
% Magnitudes of the fundamental component and harmonics
V1=100;
V3=8;
V5=15;
V7=5;

I1=100;
I3=20;
I5=15;
I7=10;
% Angles of the fundamental component and harmonics
phiv1=deg2rad(0);
phiv3=deg2rad(-70);
phiv5=deg2rad(140);
phiv7=deg2rad(20);

phii1=deg2rad(-30);
phii3=deg2rad(-165);
phii5=deg2rad(233);
phii7=deg2rad(288);
% instantaneous values
v1=sqrt(2)* V1 *sin(  w*t + phiv1);
v3=sqrt(2)* V3 *sin(3*w*t + phiv3);
v5=sqrt(2)* V5 *sin(5*w*t + phiv5);
v7=sqrt(2)* V7 *sin(7*w*t + phiv7);
v=v1+v3+v5+v7;

i1=sqrt(2)* I1 *sin(  w*t + phii1);
i3=sqrt(2)* I3 *sin(3*w*t + phii3);
i5=sqrt(2)* I5 *sin(5*w*t + phii5);
i7=sqrt(2)* I7 *sin(7*w*t + phii7);
i=i1+i3+i5+i7;

% theoretical RMS, calculated from magnitudes
V=sqrt(V1^2+V3^2+V5^2+V7^2);
I=sqrt(I1^2+I3^2+I5^2+I7^2);
% actual RMS, calculated from signals
V_RMS=sqrt(1/length(t)*sum(v.*v));
I_RMS=sqrt(1/length(t)*sum(i.*i));


P1=1/length(t)*sum(v1.*i1); % Fundamental active power (W) 3.1.2.4
P3=1/length(t)*sum(v3.*i3);
P5=1/length(t)*sum(v5.*i5);
P7=1/length(t)*sum(v7.*i7);

Q1=V1 * I1 *sin(phiv1-phii1); % Fundamental reactive power (var) 3.1.2.6
Q3=V3 * I3 *sin(phiv3-phii3);
Q5=V5 * I5 *sin(phiv5-phii5);
Q7=V7 * I7 *sin(phiv7-phii7);
%  Cross distortion powers
D13=V1 * I3;
D15=V1 * I5;
D17=V1 * I7;
D31=V3 * I1;
D51=V5 * I1;
D71=V7 * I1;
%
D35=V3 * I5;
D37=V3 * I7;
D53=V5 * I3;
D57=V5 * I7;
D73=V7 * I3;
D75=V7 * I5;
%
P=P1+P3+P5+P7;
PH=P-P1; % Harmonic active power (nonfundamental active power) (W) 3.1.2.5
QB=Q1+Q3+Q5+Q7; % Reactive power (Budenau definition) (not used in 1459)
THDV=sqrt( (V/V1)^2-1 ); % Total voltage harmonic distortion (THDV) 3.1.2.1
THDI=sqrt( (I/I1)^2-1 ); % Total current harmonic distortion (THDI) 3.1.2.1
S=V*I; % Apparent power (VA) 3.1.2.7
S1=V1*I1; % Fundamental apparent power (VA) 3.1.2.8
SN=sqrt(S^2-S1^2); % Nonfundamental apparent power (VA) 3.1.2.9
DI=S1*THDI; % Current distortion power (var) 3.1.2.10
DV=S1*THDV; % Voltage distortion power (var) 3.1.2.11
SH=S1*THDI*THDV; % Harmonic apparent power (VA) 3.1.2.12
DH=sqrt(SH^2-PH^2); % Harmonic distortion power (var)
PF1=P1/S1; % Fundamental power factor 3.1.2.15
PF=P/S; % Power factor 3.1.2.16

% draw figures
figure;
plot(t,v,'b');
hold on;
plot(t,i,'r');
grid;
xlim([0 0.1]);
xlabel('Time, s');
ylabel('Value');
legend('v','i');
% draw phasor diagram
V1_phasor=V1*exp(1i*phiv1);
V3_phasor=V3*exp(1i*phiv3);
V5_phasor=V5*exp(1i*phiv5);
V7_phasor=V7*exp(1i*phiv7);

I1_phasor=I1*exp(1i*phii1);
I3_phasor=I3*exp(1i*phii3);
I5_phasor=I5*exp(1i*phii5);
I7_phasor=I7*exp(1i*phii7);

figure;
compass(V1_phasor,'b'); hold on;
compass(V3_phasor,'g');
compass(V5_phasor,'r');
compass(V7_phasor,'k');

compass(I1_phasor,'b--');
compass(I3_phasor,'g--');
compass(I5_phasor,'r--');
compass(I7_phasor,'k--');
legend('V1','V3','V5','V7','I1','I3','I5','I7');

% draw apparent power and its components tree
figure;
axis([0 1 0 1]);
axis off;

text(0,1,['S=' num2str(S,'%6.2f') ' VA']);
annotation('arrow',[0.2 0.3],[0.9 0.8]);
annotation('arrow',[0.195 0.3],[0.9 0.68]);
text(0.23,0.83,['S_1=' num2str(S1,'%6.2f') ' VA']);
text(0.23,0.69,['S_N=' num2str(SN,'%6.2f') ' VA']);

annotation('arrow',[0.49 0.55],[0.8 0.8]);
annotation('arrow',[0.49 0.55],[0.81 0.88]);

text(0.55,0.93,['P_1=' num2str(P1,'%6.2f') ' W']);
text(0.55,0.83,['Q_1=' num2str(Q1,'%6.2f') ' var']);

text(0.55,0.69,['D_I=' num2str(DI,'%6.2f') ' var']);
text(0.55,0.59,['D_V=' num2str(DV,'%6.2f') ' var']);
text(0.55,0.49,['S_H=' num2str(SH,'%6.2f') ' VA']);

text(0.85,0.49,['P_H=' num2str(PH,'%6.2f') ' W']);
text(0.85,0.39,['D_H=' num2str(DH,'%6.2f') ' var']);

annotation('arrow',[0.49 0.55],[0.69 0.69]);
annotation('arrow',[0.49 0.55],[0.68 0.60]);
annotation('arrow',[0.49 0.55],[0.67 0.52]);

annotation('arrow',[0.72 0.78],[0.52 0.52]);
annotation('arrow',[0.72 0.78],[0.51 0.43]);

How to Evaluate Voltage Unbalance Using Symmetrical Components Method, and Determine Magnitudes by Known K2U (K0U)

To evaluate voltage unbalance we can use Symmetrical components method proposed by Charles Fortescue. Vector of three phase voltages can be expressed as a sum of three vectors (phasors): zero- positive- and negative-sequence.

\[ U_{abc} = \begin{bmatrix} U_a \ U_b \ U_c \end{bmatrix} = \begin{bmatrix} U_{a,0} \ U_{b,0} \ U_{c,0} \end{bmatrix} + \begin{bmatrix} U_{a,1} \ U_{b,1} \ U_{c,1} \end{bmatrix} + \begin{bmatrix} U_{a,2} \ U_{b,2} \ U_{c,2} \end{bmatrix} \] \[ \alpha \equiv e^{\frac{2}{3}\pi i} \] \[ \begin{align} U_{abc} &= \begin{bmatrix} U_0 \ U_0 \ U_0 \end{bmatrix} + \begin{bmatrix} U_1 \ \alpha^2 U_1 \ \alpha U_1 \end{bmatrix} + \begin{bmatrix} U_2 \ \alpha U_2 \ \alpha^2 U_2 \end{bmatrix} = \ \end{align} \] \[ \begin{align} &= \textbf{A} U_{012} \end{align} \]

From these phasors we can derive negative-sequence voltage unbalance ratio K2U and zero-sequence voltage unbalance ratio K0U

\[ K_{2U} = \frac{U_2}{U_1} \] \[ K_{0U} = \frac{U_0}{U_1} \]

Here is MATLAB script that produces the desired level of both negative- and -zero voltage unbalance ratio, then sets respective phasor magnitudes. After that we estimate unbalance using Symmetrical components method and draw a phasor diagramm.

K2U and K0U calculation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
clear all; clc;
% K2U and K0U calculator
% in three phase system we accept that UB lags UA, and UC leads UA
% in balanced three phase system all angles are 120 deg=2*pi/3 rad
% color coding: UA=red UB=green UC=blue
K2U_ref=0.05; % set K2U = 5%
% magnitudes of the three phase voltages (rms)
mag_UA=100;
mag_UB=mag_UA*(1+2*K2U_ref)/(1-K2U_ref); % set magnitude so that K2U=K0U=K2U_ref
mag_UC=100;
% angles of voltages
phi_A  = degtorad(0);    % absolute angle of UA phasor
phi_AB = degtorad(-120);  % relative angle between UA phasor and UB phasor
                         % = phi_B-phiA 
phi_AC = degtorad(120);   % relative angle between UA phasor and UC phasor
                         % = phi_C-phiA

UA=mag_UA * exp(1i*(phi_A));
UB=mag_UB * exp(1i*(phi_A+phi_AB));
UC=mag_UC * exp(1i*(phi_A+phi_AC));

% calculate symmetrical components
a=exp(1i*2*pi/3); % complex angle constant
a2=a^2;
PSP = 1/3* (1*UA  + a *UB  + a2*UC); % Positive sequence phasor
NSP = 1/3* (1*UA  + a2*UB  + a *UC); % Negative sequence phasor
ZSP = 1/3* (1*UA  + 1 *UB  + 1 *UC); % Zero sequence phasor
% calculate K2U and K0U based on symmetrical components
K2U=abs(NSP)/abs(PSP)*100 % negative sequence ratio
K0U=abs(ZSP)/abs(PSP)*100 % zero sequence ratio

ROTATE_UI=angle(UA); % draw phasors relative to UA
figure;
h=compass(UA*exp(-1i*ROTATE_UI),'r'); hold on;
set(h, 'LineWidth',2);
h=compass(UB*exp(-1i*ROTATE_UI),'g');
set(h, 'LineWidth',2);
h=compass(UC*exp(-1i*ROTATE_UI),'b');
set(h, 'LineWidth',2);
legend('UA','UB','UC');

To generate signal with known level of unbalance $K_{2U}$ we can change magnitude of one of the voltages in the balanced three phase system. So we set all angles between voltages = 120 deg. Magnitudes of Ua and Uc should be equal (Ua=Uc). Ub should be set as

\[ U_b = U_a \frac{1+2 K_{2U} }{1-K_{2U}} \]

This equation can be easily derived if we notice that

\[ U_1 K_{2U} = \frac{1}{3} (U_b + 2 U_a) K_{2U} = \frac{1}{3} (U_b - U_a) = U_2 \]

Drawing Horizontal and Vertical Lines in MATLAB Plot

Often we need to draw a vertical or horizontal line across the plot, to mark the boundaries or min/max value of the function. It’s not obvious at the first glance how to do it, becase the X and Y limits can be different, and also we need to handle zoom commands. I’m now using the solution, based on graph2d.constantline which was found here

Drawing horizontal and vertical lines in MATLAB plot
1
2
3
4
5
6
7
   GRAPH_COLOR=[.7 .7 .7];
   hx = graph2d.constantline(t_start, 'LineStyle',':', 'LineWidth',2, 'Color',GRAPH_COLOR);
   changedependvar(hx,'x');
   hx = graph2d.constantline(t_start+duration, 'LineStyle',':','LineWidth',2, 'Color',GRAPH_COLOR);
   changedependvar(hx,'x');
   hy = graph2d.constantline(magnitude, 'LineStyle',':','LineWidth',2, 'Color',GRAPH_COLOR);
   changedependvar(hy,'y');

It looks like this, and it also correctly handles zoom in and out commands.

UPD 03.04.2015 Unfortunately this does not work in Matlab 2014b or higher. Try this solution .

Fast CSV Values Export to File. MEX-based Dlmwrite Replacement for MATLAB

MATLAB functions that export values into a .CSV file are very slow. It is said, that their inefficiency is rooted in bad implementation of the string concatenation. It is refered as Schlemiel the Painter’s algorithm.

Nevertheless, sometimes I have to write huge amounts of text data into CSV files. For example I had to generate test signals for power quality measurements algoritms, such as Flicker estimation. So I had written a MEX function to speed up this task. Write process can also be interrupted by pressing Ctrl-C using this wonderful trick (http://www.caam.rice.edu/~wy1/links/mex_ctrl_c_trick/).

Function takes four parameters

mex_WriteMatrix
1
2
3
4
5
6
7
8
9
10
Usage:
     mex_WriteMatrix(filename,matrix,format,delimiter, writemode);

Parameters:
     filename  - full path for CSV file to export
     matrix    - matrix of type 'double' values to be exported
     format    - format of export (sprintf) , e.g. '%10.6f'
     delimiter - delimiter, for example can be ',' or ';'
     writemode - write mode 'w+' for rewriting file 'a+' for appending

You can find the code and compiled mexw64 for function mex_WriteMatrix here

Comparing functions (150 Mb CSV file)
1
2
 mex_WriteMatrix(fname,Z,'%10.10f',',','w+');
 dlmwrite(fname, Z, 'delimiter', ',','precision', '%10.10f');

As you can notice difference in speed is huge! On big files ~ 1GB in size, write speed values are about 17-18 Megabytes/s.

You can use it freely, I hope it will save you a lot of time.

Using Coursera-dl to Download Lecture Videos and Materials From Coursera.org

I use Coursera-dl to download lecture videos and materials from Coursera, so I can watch them later after the course is finished. How to install coursera-dl on Windows machine? It’s easy.

1.First of all you need to download Python 2.7 (don’t download Python 3, it is a different and incompatible version) from Python download page.

2.Then install Pip (Python package manager)

  • Download file https://bootstrap.pypa.io/get-pip.py to Python directory.
  • Run the command line cmd
  • Change directory to your Python installation ( C:\Python27 for my case )
  • execute python get-pip.py

3.Install Coursera-dl script

  • In the command line change directory to pip (C:\Python27\Scripts for my case)
  • execute pip install coursera-dl

4.Now you are ready to download materials! For this you need the class name identificator from the browser URL.

Provide your user name and password, directory to download and class name identificator as the parameters into coursera-dl script. Don’t forget, that you have to accept the honor code of the class before using this script (happens the very first time you go to the class page).

Sample command line to download Coursera materials for current Machine Learning class
1
coursera-dl -u myusername -p mypassword -d C:\Coursera\ ml-006m