Применение пакетного вейвлет-преобразования для определения составляющих мощности при несинусоидальных режимах Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

Оригинальная статья / Original article УДК: 621.316

DOI: 10.21285/1814-3520-2016-8-136-145

ПРИМЕНЕНИЕ ПАКЕТНОГО ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ ДЛЯ ОПРЕДЕЛЕНИЯ СОСТАВЛЯЮЩИХ МОЩНОСТИ ПРИ НЕСИНУСОИДАЛЬНЫХ РЕЖИМАХ

© Л.А. Файфер1, Д.С. Осипов2, Е.Н. Еремин3, Н.Н. Долгих4

Омский государственный технический университет, 644043, Россия, г. Омск, пр-т Мира, 11.

Резюме. Цель. Для исследования показателей качества электроэнергии при несинусоидальных режимах большое значение имеет правильное определение их параметров: действующих значений токов, напряжений, активной и реактивной мощности. Целью данной работы является модернизация существующих алгоритмов расчета мощностей при несинусоидальных режимах с применением пакетного вейвлет-преобразования. Методы. Показан способ определения действующего значения тока несинусоидального режима с предварительной частотной декомпозицией исследуемого сигнала с использованием пакетного вейвлет-преобразования - современного метода цифровой обработки сигнала, позволяющего сохранять информацию о сигнале не только в амплитудно -частотном, но и во временном пространстве. Результаты. Модернизированы существующие математические модели для расчета мгновенной, активной, реактивной мощности и мощности искажения при несинусоидальных режимах на основе пакетного вейвлет-преобразования. Проведено сравнение методики определения параметров по Фризе и аналитического расчета. Заключение. Решена задача частотной декомпозиции параметров режима электроэнергетической системы. Подход на основе пакетного вейвлет-разложения позволяет сужать частотные диапазоны, а предлагаемая методика повышает точность расчетов несинусоидальных режимов при непрерывном контроле.

Ключевые слова: пакетное вейвлет-преобразование, несинусоидальный режим, активная мощность, реактивная мощность, теория Фризе, преобразование Фурье.

Формат цитирования: Файфер Л.А., Осипов Д.С., Еремин Е.Н., Долгих Н.Н. Применение пакетного вейвлет-преобразования для определения составляющих мощности при несинусоидальных режимах // Вестник Иркутского государственного технического университета. 2016. № 8 (115). С. 136-145. DOI: 10.21285/1814-3520-2016-8-136-145

IMPLEMENTATION OF WAVELET PACKET TRANSFORM FOR POWER COMPONENTS DETERMINATION

IN NON-SINUSOIDAL SITUATIONS

L.A. Faifer, D.S. Osipov, E.N. Eremin, N.N. Dolgikh

Omsk State Technical University, 11 Mira pr., Omsk, 644043, Russia.

Abstract. Purpose. The study of power quality indices in non-sinusoidal modes is largely dependent on correct determination of their parameters: actual values of currents, voltages, active and reactive power. The purpose of this work is to upgrade the existing algorithms of power calculation in non-sinusoidal situations with the application of wavelet packet transform. Methods. The paper shows a method for determining the actual value of current in the non-sinusoidal mode with the frequency predecomposition of the tested signal and the use of the wavelet packet transform. The latter is a modern method of digital processing of a signal allowing to store information about the signal both in the amplitude and frequency space and in the time space. Results. Based on the wavelet packet transform the existing mathematical models for the calculation of instantaneous, active, reactive power and distortion power in non-sinusoidal modes have been upgraded. The Fryze method of parameter estimation has been compared with the analytical calculation. Conclusion. The problem of frequency decomposition of power system mode parameters is solved. A wavelet packet decomposition-

1

Файфер Лилия Андреевна, магистрант, e-mail: faiferlilia@mail.ru Faifer Liliya, Master's degree student, e-mail: faiferlilia@mail.ru

2Осипов Дмитрий Сергеевич, кандидат технических наук, доцент кафедры электроснабжения промышленных предприятий, e-mail: ossipovdmitriy@list.ru

Osipov Dmitriy, Candidate of technical sciences, Associate Professor of the Department of Power Supply of Industrial Enterprises, e-mail: ossipovdmitriy@list.ru

3Еремин Евгений Николаевич, доктор технических наук, профессор кафедры машиностроения и материаловедения, e-mail: weld_techn@mail.ru

Eremin Evgeniy, Doctor of technical sciences, Professor of the Department of Mechanical Engineering and Materials Science, e-mail: weld_techn@mail.ru

4Долгих Надежда Николаевна, аспирант, e-mail: nabal2006@list.ru Dolgikh Nadezhda, Postgraduate, e-mail: nabal2006@list.ru

based approach allows to narrow the frequency ranges while the proposed technique improves the calculation accuracy of non-sinusoidal modes under continuous monitoring.

Keywords: Wavelet Packet transform, non-sinusoidal situation, active power, reactive power, Fryze theory, Fourier transform

For citation: Faifer L.A., Osipov D.S., Eremin E.N., Dolgikh N.N. Implementation of wavelet packet transform for power components determination in non-sinusoidal situations. Proceedings of Irkutsk State Technical University. 2016, no. 8 (115), pp. 136-145. DOI: 10.21285/1814-3520-2016-8-136-145

Введение

Стратегия развития электросетевого комплекса Российской Федерации на фоне ежегодного увеличения спроса на электрическую энергию обусловливает необходимость внедрения новых технологий и повышения энергоэффективности энергетической отрасли. Ставится задача внедрения технологий, использующихся в сетевых компаниях развитых стран, в частности, технологии «умных» электрических сетей, позволяющих «повысить пропускную способность и стабильность энергоснабжения, сократить потери и издержки на технический и коммерческий учет у потребителя»5.

Это предопределяет актуальность развития научных исследований методов цифровой обработки и анализа режимных параметров электроэнергетических систем (ЭЭС), а также разработки методов математического моделирования установившихся и переходных процессов в условиях неопределенности исходной информации, связанной с вариациями значений вырабатываемых и потребляемых мощностей и конфигурации системы.

В настоящее время получает распространение метод цифрового анализа сигналов на основе вейвлет-преоб-

5

Стратегия развития электросетевого комплекса Российской Федерации: Постановление Правительства РФ от 3 апреля 2013 г. № 511-р // Собрание законодательства РФ. 2013. № 14. Ст. 1738. Development strategy of the national grid complex of the Russian Federation: Government resolution of the Russian Federation of April 3, 2013, no. 511 -p. Collection of the legislation of the Russian Federation. 2013, no. 14, article 1738.

6Яковлев А.Н. Введение в вейвлет-преобразование: учеб. пособие. Новосибирск: Изд-во НГТУ, 2003. 104 с.

Jakovlev A.N. Introduction to wavelet transform, learning aids. Novosibirsk, NGTU Publ., 2003, 104 p._

разования (ВП). В [1] проведен анализ зарубежных источников, определивший наиболее популярные направления использования ВП в ЭЭС: диагностика электрооборудования; анализ переходных процессов в ЭЭС; релейная защита и автоматика ЭЭС и др. Преимущества ВП по сравнению с преобразованием Фурье (ПФ) при анализе нестационарных сигналов определили перспективность ВП в исследовании несинусоидальных режимов электроэнергетических систем. Алгоритмы для непрерывного контроля параметров качества электроэнергии в режиме реального времени на основе ВП представлены в работе.

Вейвлет-преобразование исследуемой функции тока (напряжения, э.д.с. и т.д.) ¡(г) по составляющим в различных масштабах и частотных диапазонах может быть реализовано путем свертки ¡(г) и вейвлет-функции у (г). Вейвлет-функция, сдвинутая относительно времени и масштаба у к(г), определяется по формуле6

у]Л(г) = -= 11(г)у(2Ч - к) Жг, (1) л/21 -м

где к и 1 - сдвиг по времени и масштабирующий коэффициент соответственно; - нормирующий показатель.

л/27

Принято различать непрерывное ВП и дискретное ВП. В данной работе проиллюстрируем частотную декомпозицию параметров несинусоидального режима на основе разновидности дискретного ВП, называемого пакетным вейвлет-преобразованием.

Теоретический базис вейвлет-преобразования

Непрерывный (аналоговый) сигнал несинусоидального тока может быть представлен рядом Фурье, т.е. суммой гармонических колебаний:

i(t) = h + Ilmsin(at -ф ) + +I2m sin( cot - ф2 ) + I3m sin( cot + ф3 ) + ....

(2)

п = Тп¥й . (4)

Для рассматриваемого примера п = 0,02 • 400 = 8 отсчетов. Следовательно, исходный дискретный сигнал тока может быть задан вектор-строкой:

i( ts) = [^

Хл i'>

Vi]-

(5)

Цифровой сигнал определяется цифровыми отсчетами ¡(^) с заданной

частотой дискретизации ^, определяемой теоремой Котельникова:

F, > F

А d — А m '

(3)

где

F..

максимальная

частота

сигнала, Гц.

На рис. 1 приведены тестовые сигналы тока в виде непрерывной и дискретной величин. Максимальная частота сигнала ^ = 150 Гц, поэтому по выражению (3)

частота дискретизации была принята ^ = 400 Гц. Частота дискретизации сигнала

^ и время Гт, в течение которого исследуется сигнал, определяют количество мгновенных отсчетов дискретной величины:

Допустим, полученная последовательность мгновенных значений является результатом вейвлет-преобразования начального уровня, который обозначим j0 = 0. Для последующего j1 = j0 + 1 = 1 уровня разложения схема пакетного вейвлет-преобразования представлена на рис. 2. В результате операции децимации (прореживания) происходит сокращение количества элементов в каждом узле (m = 0, m = 1) в два раза. Таким образом, мы получаем 2 узла пакетного вейвлет-разложения. Исходная числовая последовательность мгновенных значений токов состояла из 8 элементов, после свертки и децимации числовая последовательность в узле j = 1, m = 0 будет состоять из 4 элементов:

ij, m(k) = ¿1, 0 (k) = ij, m(k) =

= [h,o(0) ( 1) ¿1,0(2) ¿1,0(3)1

(5)

Рис. 1. Непрерывный и дискретный сигналы тока Fig. 1. Continuous and discrete signals of current

Рис. 2. Схема пакетного вейвлет-разложения Fig. 2. Diagram of wavelet packet decomposition

В общем случае k - порядковый номер элемента числовой последовательности. Для k можно записать следующее неравенство:

Ыш < ^ < Ы(ш +1) -1, (6)

где Ы - уровень разложения; ш - номер узла пакетного вейвлет-преобразования.

Для реализации вейвлет-разложения необходимо задать масштабирующую функцию ф(1) и вейвлет-функцию у/(1), которые генерируют базисы, используемые для разложения и реконструкции сигнала. Одними из простейших вейвлет-функций являются вейвлеты Хаара, которые могут быть заданы по выражениям:

Вейвлеты Хаара, полученные по (7), (8), представлены на рис. 3.

С применением вейвлета Хаара для уровня разложения j = 1, и узлов m = 0 и m = 1 аппроксимирующие и детализирующие вейвлет-коэффициенты будут получены следующим образом:

, /m-i0°(0) + V- 0 +1 0607 'l,0 ' 0)

42

h, i( 4) = ^

42

il, 0 (k) =

1, 0 < t < 1 / 2

w(t) = J_1, 1 / 2 < t < 1

<P(t) = ■

" 42

с/2) 0 _ 1,0607

= 42

i1,0 ( 0) i1,0 ( 1) i1,0 ( 2) i1,0 (3) 0,7500 1,1036 = _0,7500 _1,1036

= 0,75;

= _0,75

(7)

0, 1 < t < 0 Ъ( 4 ) _0, 7500

1, 0 < t < 1 (8) i11(k) = iu( 5) 1( 6) = _0,3964 0,7500

0, 1 < t < 0. h,, / 7) 0,3964

-1

V2

У2

Рис. 3. Вейвлеты Хаара Fig. 3. Haar wavelets

Действующее значение несинусоидального тока, представленного на рис. 1, может быть определено через вейвлет-коэффициенты:

1 j,m

2 Jn(m+1 )-1

1 S

i.

/к) .

(9)

2 J nm

Если требуется провести частотную декомпозицию сигнала, т.е. выделить из результирующего сигнала действующее значение тока определенной полосы частот, то необходимо в формуле (9) использовать вейвлет-коэффициенты уровня разложения j и узла m, которые отвечают за искомый диапазон частот. В нашем случае для проведения численного эксперимента сигнал тока представлен суммой двух частот - 50 Гц и 150 Гц:

i(t) = 1 • sin( 2 -ж- 50 • t) +

+ 0,5 • sin(2-ж-150 • t).

Таким образом, для нахождения действующего значения тока основной частоты достаточно первого уровня (j = 1) вейвлет-разложения (рис. 2) для узла m = 0:

I( 0< f <100 ; (if/ 0) + &( 1) + ¿i,0 ( 2 ) + 3)) =

= . 11 (0,75002 +1, 0362 + (-0, 75 00f + (-1, 036/ ) =

= 0, 6395 (А).

Аналогично для уровня вейвлет-разложения } = 1 и узла т = 1 может быть получено действующее значение тока третьей гармоники:

' (100,/ <200) = 1 (¡и ( 4) + ¡121(5) + ¡и ( 6) + ¡и ( 7)) =

= 0, 4242 (А).

Погрешность вычисления действующих значений отдельных гармоник через вейвлет-коэффициенты относительно фактических значений, определяемых через амплитудное значение, обусловлена выбором типа вейвлета. Вейвлет Хаара является достаточно удобным для иллюстрации принципов ВП, но для задач электротехники демонстрирует существенную погрешность [2]. Для рассматриваемого численного примера проведем процедуру пакетного вейвлет-разложения при помощи других типов вейвлетов. Как следует из табл. 1, достаточно приемлемые результаты в определении действующего значения тока отдельных частотных компонент демонстрирует вейвлет Ингрид Добеши 42-го порядка, поэтому в дальнейшем будем использовать этот тип вейвлета.

Таблица 1

Сравнение точности вычисления действующего значения тока различными типами вейвлетов

Table 1

Comparison of calculation accuracy of actual current value by various types of wavelets

Аналитический расчет / Analytical calculation Расчет по вейвлет-коэффициентам для вейвлетов / Calculation by wavelet-coefficients for wavelets

Хаара / Haar Добеши 10-го порядка / 10th order Daubechies Добеши 42- го порядка / 42d order Daubechies

Действующее значение основной частоты /Г50), А / Actual value of the main frequency /Г50), A 0,7071 0,6395 0,7095 0,7071

Действующее значение третьей гармоники Iri50А / Actual value of the third harmonic Iri50А 0,3536 0,4242 0,3487 0,3535

Погрешность определения действующего значения, I(50) /1(150), % / Error of actual value determina- ^ I(50) / I( 150) ,% - 9,6 % / 20,0 % 3,4 % / 1,4 % менее / less than 0,1 %

Методика расчета мощности по вейвлет-коэффициентам

Обзор различных подходов применения вейвлет-преобразования для расчета составляющих мощности при несинусоидальных режимах на основе дискретного ВП представлен в [3]. Однако дискретное ВП имеет логарифмическую шкалу частотных коридоров, что позволяет с успехом решать задачи по разложению исходного сигнала на две компоненты: основную частоту и высшие гармоники. Выделить же различные полосы частот (высших гармоник) в отдельные составляющие позволяет пакетное ВП.

Целью данной работы является модернизация существующих алгоритмов расчета мощностей при несинусоидальных режимах с применением пакетного ВП. Напряжение и ток могут быть представлены коэффициентами, полученными в результате пакетного вейвлет-разложения:

2 Jn(m+1 )_1

2 Jn(m+1 )_1

i(t) = S iJ,m (k) -w(t)

(10)

u(t) = S Ujm(k) -V(k) ■

(11)

Активная мощность через вейвлет-коэффициенты для искомого частотного диапазона может быть определена как

( 2-3п(ш+1 )-1 ^

P = 1

J,m

J n

ij,m(k)Uj,m(k)

(12)

Определение реактивной мощности при несинусоидальном режиме по-прежнему является обсуждаемой проблемой. На этот счет нет однозначного определения, так как общепринятое определение реактивной мощности при наличии гармоник теряет смысл. Тем не менее, существует несколько различных подходов к решению данного вопроса: определение реактивной мощности по Буденау, по Шарону и Кастерсу-Муру, по Маевскому и др.7 Захарова Т.В. Вейвлет-анализ и его приложения: учеб. пособие. М: Инфра-М, 2012. 158 с. Zaharova T.V. Wavelet analysis and his applications: learning aids, Moscow, Infra-M Publ., 2012, 158 p.

2 nm

2

nm

В исследовании не ставится задача обоснования справедливости той или иной теории. Цель проводимого исследования -показать возможность применения пакетного вейвлет-преобразования для нахождения составляющих мощности по теориям, предложенным ранее ведущими учеными в данной области.

По теории реактивной мощности Фризе полный ток состоит из двух составляющих: активного и реактивного тока [4]:

I = I + jI .

a J r

(13)

Реактивная мощность по Фризе определяется по формуле

Q = UI, =4UI2-UI2a .

(14)

В соответствии с поставленной задачей адаптируем теорию реактивной мощности по Фризе для алгоритма на основе пакетного ВП. Для выделения из полного тока активной составляющей необходимо предварительно вычислить активную проводимость исследуемой ветви по формуле

P,

g =

(Um)2

1

( a" Jn(m+1 )-1

n

_ V

S ij. m(k)uj m m(k)

a j nm

(13)

1 a J n(m+1 )-1

1 S ujm(k) n J

Активную составляющую мгновенного тока вычислим по формуле

I = ей. =

а о т

( 2-]п(т+1)-1 Л

I 1Ьт(к)иьт(к)

1

n

__V

a- Jr

1 a J n(m+1 )-1

1 S ujm(k)

n a- J nm 1 a- J n( m+1 )-1

- S U,m(k).

(14)

a- jn

Действующее значение активного тока определяем как

I =

1 a-J n( m+1 )-1

TW S lKk)

a J nm

1

1 a - J n( m+1 )-1

1 S i2 (k).

■>W ¿J J. m ' /

(15)

a- Jn

Результаты численного моделирования

Дли иллюстрации эффективности предлагаемой методики расчета параметров режима по вейвлет-коэффициентам пакетного вейвлет-преобразования прове-

дем численный эксперимент для системы электроснабжения, представленной на рис. 4.

Система

ТМ 40/10 СШ 0,4кВ

СП

Параметры элементов расчетной схемы / Parameters of computational scheme elements

Трансформатор ТМ 40/10 / Transformer TM 40/10 S = 40 кВА, Ub = 10 кВ, Uh = 0,4 кВ, Rt = 88 мОм, Хт =157 мОм

Кабель (КЛ) АВВГ 3x50+1x35 / Cable (CL) AVVG 3X50+1x35 1дл.доп = 110 А Rü = 0,625 мОм/м, Хй = 0,085 мОм/м, L = 200 м

Нагрузка (Н) / Loading (L) S = 20 кВА, cos ф = 0,95

Выпрямитель (В) / Rectifier (R) S = 13,7 кВА

Рис. 4. Расчетная схема Fig. 4. Computational scheme

х

Н

Таблица 2

Результаты расчета параметров несинусоидального режима

Table 2

Results of non-sinusoidal mode parameters calculation_

Параметр/ Parameter Частота, Гц / Frequency, HZ Обозначение, единица измерения / Designation, unit of measure Мгновенное значение / Instantaneous value

Напряжение на СШ 0,4 кВ / Voltage on СШ of 0.4 kV 50 u, В 220V2sin(or)

150 u3, В 4V2sin(5^r)

250 щ, В 3V2sin(7^t)

Ток в кабельной линии / Current in a cable line 50 h, В 30V2sin(>r -18°)

150 h, в 5V2sin(5^r - 44°)

250 ¿5 , В 3V2sin(7^r - 60°)

150 200 250

Частота (Гц) Frequency (Hz)

Рис. 5. Амплитудно-частотное представление напряжения на СШ 0,4 кВ Fig. 5. Amplitude-frequency representation of voltage on СШ of 0.4 kV

so loo 150 гоо 250 зоо 350 joo Частота (Гц) Frequency (Hz)

Рис. 6. Амплитудно-частотное представление тока в кабеле Fig. 6. Amplitude-frequency representation of current in a cable

По представленной выше методике произведем расчет параметров режима с разложением сигналов тока и напряжения пакетным вейвлет-преобразованием. Частотные диапазоны вейвлет-разложения соответствуют схеме, представленной на

рис. 2. Итоговые результаты, а также сравнение с аналитическим расчетом параметров режима представлены в табл. 3. Истинные значения вычислены через действующие значения токов и напряжений.

Расчет параметров режима по вейвлет-коэффициентам Calculation of mode parameters by wavelet coefficients

Таблица 3 Table 3

Параметр/ Parameter Частота, Гц / Frequency, Hz Обозначение, единица измерения / Designation, unit of measure Истинные значения / True values Значение, вычисленное через вейвлет-коэффициенты / Value, calculated by wavelet coefficients Погрешность, % / Error, %

Активная составляющая тока / Active component of current 50 ^, А 28,53169 28,5317 0,00003

150 , А 3,59669 3,5967 0,0003

250 ^, А 1,5 1,5 0

Активная мощность / Active power 50 Р, Вт 6277,9 6277,8 0,0016

150 Р, Вт 14,387 14,386 0,0048

250 Р, Вт 4,5 4,353 3,2578

Реактивная мощность и мощность искажения / Reactive power and distortion power 50 + N, вар 2039,5 2032,1 0,3643

150 Q + N, вар 13,8932 13,8765 0,1203

250 Q + N, вар 7,7942 7,4905 4,0545

Заключение

На основе рассмотренной методики путем применения пакетного вейвлет-преобразования решена задача частотной декомпозиции параметров режима электроэнергетической системы. По сравнению с применяемым для решения аналогичных задач дискретным вейвлет-преобра-зованием, предлагаемый подход значи-

тельно сужает частотные диапазоны, что позволяет с достаточной степенью точности идентифицировать и определять токи, напряжения и мощности для каждой из частот высших гармоник при несинусоидальных режимах. Кроме того, повышается точность расчета несинусоидальных режимов при непрерывном контроле.

Библиографический список

1. Мисриханов А.М. Применение методов вейвлет-преобразования в электроэнергетике // Автоматика и телемеханика. 2006. № 5. С. 5-23.

2. Долгих Н.Н. Определение коэффициентов искажения синусоидальности формы кривой тока по вейвлет-коэффициентам // Политематический сете-

вой электронный научный журнал Кубанского государственного аграрного университета. 2015. № 113. С. 814-825.

3. Morsi W.G. On the implementation of time-frequency transforms for defining power components in non-

sinusoidal situations. Electric Power Components and Systems, 2009, vol. 37, no. 4, pp. 373-392. 4. Жемеров Г.Г., Ильина О.В. Теория мощности Фризе и современная теория мощности // Електро-технка I Електромеханка. 2007. № 6. С. 63-65.

References

1. Misrihanov A.M. Primenenie metodov vejvlet-preobrazovanija v jelektrojenergetike [Application of wavelet transform methods in electric power industry]. Avtomatika i telemehanika [Automation and Remote Control]. 2006, no. 5, pp. 5-23. (In Russian)

2. Dolgih N.N. Opredelenie kojefficientov iskazhenija sinusoidal'nosti formy krivoj toka po vejvlet kojefficien-tam [Determining total harmonic distortion of current by wavelet coefficients]. Politematicheskij setevoj jel-ektronnyj nauchnyj zhurnal Kubanskogo gosudarstven-nogo agrarnogo universiteta [Polythematic on-line elec-

Критерии авторства

Файфер Л.А., Осипов Д.С., Ерёмин Е.Н., Долгих Н.Н. имеют равные авторские права; Долгих Н.Н. несет ответственность за плагиат.

Конфликт интересов

Авторы заявляют об отсутствии конфликта интересов.

Статья поступила 02.06.2016 г.

tronic scientific journal of the Kuban state agricultural university]. 2015, no. 113, pp. 814-825. (In Russian)

3. Morsi W.G. On the implementation of time-frequency transforms for defining power components in non-sinusoidal situations. Electric Power Components and Systems, 2009, vol. 37, no. 4, pp. 373-392. (In English)

4. Zhemerov G.G., Il'ina O.V. Teorija moshnosti Frize i sovremennaja teorija moshhnosti [Fryze theory of power and modern theory of power]. Elektrotehnika I El-ektromehanika [Electrical engineering and electrome-chanics]. 2007, no. 6, pp. 63-65. (In Russian)

Authorship criteria

Faifer L.A., Osipov D.S., Eremin E.N., Dolgikh N.N. have equal author's rights. Dolgikh bears the responsibility for plagiarism.

Conflict of interests

The authors declare that there is no conflict of interest regarding the publication of this article.

The article was received on 02 June 2016