Основное
Передаточная функция резонансного регулятора в непрерывной форме. Коэффициент Kr имеет тот же смысл что интегральная составляющая в ПИ-регуляторе
Резонансный регулятор. Kr играет роль коэффициента Ki интегрального регулятора для постоянных сигналов Резонансный регулятор о��еспечивает бесконечно большой коэффициент усиления на частоте f₀, обеспечивая астатическое регулирование одновременно по амплитуде и фазе сигнала с этой фиксированной частотой.
Что это такое и для чего необходимо
Иногда возникает необходимость стабилизировать напряжение или ток преобразователя, например, бесперебойного источника питания, солнечного инвертора и иных устройств где на выходе имеется синусоида.
Для трёхфазных систем обычно используется преобразование Парка-Горева. Вначале трёхфазный сигнал ABC преобразуется в ортогональный αβ0 (плюс нулевая последовательность) посредством преобразования Кларке, который, затем, посредством преобразованием Парка‑Горева, становится проекциями на синхронно вращающуюся систему координат, представляющие собой сигналы с постоянной составляющей, содержащие информацию о фазе и амплитуде исходной трёхфазной системы. В идеальном варианте симметричной системы, содержащей только прямую последовательность, эти сигналы не содержат вторую гармонику. Управление по симметричным компонентам — отдельная тема, следует отметить, что резонансный регулятор позволяет реализовать независимое управление по отдельным фазам.
В однофазных системах используется восстановление ортогональной составляющей производной/интегратором/фазовым фильтром либо специальные методы по измерению амплитуды и относительной фазы требуемых величин применяя счётчики и компараторы, однако это не всегда надёжно и «плавно». В этом случае на помощь приходит резонансный регулятор, где изначально в качестве сигнала задания используется синусоидальный сигнал с необходимой амплитудой и фазой.
Анализ замкнутой системы с Р-регулятором в непрерывной форме
Общая структура представлена на рис. 1. Имеется входной синусоидальный сигнал, ПР-регулятор, объект управления в общем виде.
С использованием предложенной методики покажем равенство фазы и амплитуды в установившемся режиме. Общий скрипт для вычислений в системе символьных вычислений Maxima:
/*Передаточная функция ПР-регулятора*/ V1:Kr*p/(p^2+(2*%pi*f_0)^2); /*Уравнение замкнутой системы*/ V2:[e=x-y,y=Fk(p)*V1*e];
Вначале задаётся вид регулятора, составляется система уравнений для замкнутой системы на рис. 1.
Далее формируется решение, подставляется в него значение входного сигнала x с целью получения переходной характеристики:
/*Решение системы относительно выходной величины*/ V3:factor(y),solve(V2,[e,y]); /*Задание входного сигнала во временной форме*/ V4:A*sin(2*%pi*f_0*t+ψ_0); /*Преобразование Лапласа от входного сигнала*/ V5:x=laplace(V4,t,p); /*Подстановка входного сигнала в решение системы уравнений - переходная характеристика*/ V6:factor(V3),V5;
Производится умножение на коэффициент перехода к комплексно-частотной точке в установившемся режиме
/*Задание коэффициента перехода к комплексно-частотной точке*/ V7:(p^2+4*%pi^2*f_0^2)/(2*%pi*f_0); /*Умножение на коэффициент и упрощение*/; V8:factor(V6*V7); /*Переход к комплексно-частотной характеристике в заданной точке*/ V9:p = %i*2*%pi*f_0; /*Комлпексное значение в установившемся режиме*/ V10:factor(V8),V9;
Производится подстановка для перехода к комплексно-частотной точке в установившемся режиме
Легко заметить что получается ничто иное как формула Эйлера. Таким образом, амплитуда выходного сигнала в установившемся режиме равна A и фаза равна ψ₀, что и является непосредственно входным сигналом. Обратить внимание что выходной сигнал не зависит от передаточной функции объекта управления в установившемся режиме по интересующей гармонике.
Чтобы это показать наглядно, зададим необходимую функцию для объекта управления. Пусть это будет апериодическое звено первого порядка вида ФНЧ с постоянной времени τ=0.1 секунды, определяется как 1/(τ•p+1).
Остальные параметры соответственно равны Kr = 1000, ψ₀ = -0.2, f₀ = 50, A = 1.3.
Далее приведены результаты численного анализа для схемы, составленной в среде Scilab+xCos.
Параметры моделирования и исходные данные
Результаты моделирования. В разомкнутой системе амплитуда сигнала достаточно мала и соответствует прохождению через ФНЧ, при этом фаза является отстающей и примерно равна 90°.
Получим тот же результат для проверки с использованием Maxima, продолжим скрипт с использованием вспомогательной библиотеки для численных методов, позволяющих осуществить обратное преобразование Лапласа для функций высших порядков.
/*Задать форму объекта управления Fk(p)*/ V11:Fk(p)=1/(τ*p+1); /*Подставить объект управления в переходную характеристику*/ V12:factor(V6),V11; /*Задать численные параметры*/ Num:[A=1.3,Kr=2000,ψ_0=-0.2,f_0=50,τ=0.02]; /*переходная характеристика в численном виде*/ V13:V12,Num; /*Загрузить пакет для работы с обратным преобразованием Лапласа в численном виде*/ load("coma"); /*получить обратное преобразование Лапласа */ V14:nilt(V13,p,t); /*Входной сигнал с численными значениями*/ V15:V4,Num; /*Построение графика переходного процесса*/ wxplot2d([V15,V14],[t,0,0.4],[color,green,black],[legend,false]),wxplot_size=[480,320];
Задание объекта управления и получение переходной характеристики
После чего можно подставить численные параметры, загрузить пакет «coma» для работы с функциями высших порядков (более 4-го). Команда по-умолчанию ilt работает до 4 порядка включительно, для преодоления этого ограничения была разработана nilt (numerical inverse Laplace transform).
Результат представлен далее, видно, что выходной сигнал стремится к необходимой фазе и амплитуде 1.3.
Таким образом, резонансный регулятор позволяет реализовать астатическое (с нулевой ошибкой) управление как по амплитуде, так и по фазе сигнала управляющего воздействия, обеспечивая в установившемся режиме их равенство.
Построим частотную характеристику регулятора, чтобы показать бесконечную добротность на частоте управления.
/*Переход к комплексно-частотной характеристике (КЧХ) для всех частот*/ V16:p = %i*2*%pi*f; /*КЧХ регулятора*/ V17:V1,V16; /*АЧХ и ФЧХ регулятора, предположить что параметры неотрицательные для упрощения не содержащего знак модуля*/; assume(Kr>0,f>0,f_0>0); /*Действительная и мнимая части КЧХ*/ V17R:realpart(V17); V17I:factor(imagpart(V17)); /*АЧХ и ФЧХ регулятора, предположить что параметры неотрицательные для упрощения не содержащего знак модуля*/; assume(Kr>0,f>0,f_0>0); V17A:factor(cabs(V17)); V17P:carg(V17); /*Подстановка численных параметров */ V18A:V17A,Num;V18P:V17P,Num; /*График АЧХ*/ wxplot2d(V18A,[f,0,400],[y,0,10],[grid2d,true],[color,green,black],[legend,false],[ylabel,"|f(j2πf)|"]),wxplot_size=[480,320]; /*Построение графика ФЧХ*/ wxplot2d(V18P*180/%pi,[f,0,400],[y,0,360],[grid2d,true],[ylabel,"∠f(j2πf)"],[color,green,black],[legend,false]),wxplot_size=[480,320]; /*Построение графика мнимой и действительной частей*/ wxplot2d([V18R,V18I],[f,0,400],[y,-10,10],[grid2d,true],[ylabel,"Re,Im f(j2πf)"],[color,green,black],[legend,false]),wxplot_size=[480,320];
Осуществляется подстановка для оператора Лапласа, осуществляющая переход в комплексно-частотную область для всех значений частоты, далее осуществляется взятие модуля, угла и нахождение действительной и мнимой частей КЧХ, обратить внимание, что КЧХ — чисто мнимая.
На рис. 14 представлены результирующие характеристики. Обратить внимание на нулевую действительную часть КЧХ и асимптоты, соответствующие частоте 50 Гц.
Цифровая форма реализации
Для получения регулятора в дискретном виде воспользуемся методом согласования нулей и полюсов. Об этом методе будет отдельная публикация равно как и о реализации в Maxima цифровых фильтров.
Для начала используем следующий код в новом скрипте
/*Инициализация*/ kill(all);ratprint:false;fpprintprec:4; /*Резонансный регулятор*/ V1:Kr*p/(p^2+(2*%pi*f_0)^2); /*Знаменатель и числитель*/ V1n:num(V1);V1d:denom(V1);
Общая схема (рис. 15) метода нулей и полюсов представлена на рисунке, здесь и далее fd - частота дискретизации
Далее находятся корни числителя и знаменателя, записывается передаточная функция в общем виде исходя из значения корней. Обратить внимание на использование функции demoivre, которая преобразует комплексные экспоненты в синусы-косинусы согласно формулы де Муавра.
/*Корни знаменателя и числителя*/ Kn:solve(V1n,p);Kd:solve(V1d,p); /*Нахождение произведения разностей вида z-exp(p[i]/fd) для i-го корня*/ /*Для числителя*/ Zn:product(subst(Kn[i],(z-exp(p/fd))),i,1,length(Kn)); /*Для знаменателя*/ Zd:trigsimp(demoivre(product(subst(Kd[i],(z-exp(p/fd))),i,1,length(Kd)))); /*Передаточная функция в общем виде*/ V2:Krd*Zn/Zd; /*Сравнение частотных характеристик*/; /*Переход в КЧХ для непрерывной функции*/ V3:p=%i*2*%pi*f; /*Переход в КЧХ для дискретной функции*/ V4:z=exp(%i*2*%pi*f/fd); /*Подстановка в непрерывную ПФ*/ V5:V1,V3; /*Подстановка в дискретную ПФ*/ V6:factor(V2),V4;
Вывод представлен на рисунке 16 в виде последовательности выражений
Следует отметить, что билинейное преобразование (метод Тастина) смещает частоты, поэтому используется именно метод согласования нулей и полюсов, дающий точные частотные точки, однако, при этом, искажается амплитуда. Для этого необходимо найти вспомогательный коэффициент Krd, который приближает непрерывную функцию к дискретной в заданной частотной точке по амплитуде.
/*Подстановка в непрерывную ПФ*/ V5:V1,V3; /*Подстановка в дискретную ПФ*/ V6:factor(V2),V4; /*Потребуем, чтобы скорость изменения составляющих КЧХ была одинаковой на нулевой частоте. Используется производная*/ DV5:subst(f=0,diff(V5,f)); DV6:subst(f=0,diff(V6,f)); /*Решение для коэффициента, удовлетворяющего равенству производных на нулевой частоте*/ V7:solve(DV5=DV6,Krd); /*Зададим численные параметры*/ Num:[Kr=2000,f_0=50,fd=1000]; /*Нахождение вспомогательного коэффициента*/ V7N:float(V7),Num; /*Общий вектор численных параметров*/ Num:append(V7N,Num); /*Упрощение выражений для АЧХ резонансного регулятора в непрерывном и дискретном виде*/ assume(f>0,f_0>0,fd>0,Kr>0,Krd>0); V5A:factor(cabs(V5));V6A:factor(trigreduce(trigsimp(factor(cabs(V6))))); /*Подстановка численных параметров*/ V5N:float(V5A),Num;V6N:float(V6A),Num; /*Построение АЧХ*/ wxplot2d([V5N,V6N],[f,0,2000],[y,0,5],[grid2d,true],[ylabel,"|Fd(p)|, |Fz(p)|"],[color,green,black],[legend,false]),wxplot_size=[480,320];
При обычном методе где на нулевой частоте отличное от нуля значение, например, для ФНЧ, вспомогательный коэффициент находится подстановкой p=0 и z=1. Если же необходимо обеспечить одинаковый рост частотной характеристики вблизи этой точки, особенно если она равна нулю, используется производная. Порядок нахождения представлен на рис. 17.
После подстановки этого коэффициента в дискретную передаточную функцию, подстановки численных параметров, можно построить следующий график АЧХ (рис. 18)
Далее реализуем цифровой фильтр и сравним результат во временной форме. Получим сперва результат обратного преобразования Лапласа от воздействия на резонансный регулятор ступенчатой функции, он довольно прост: (20•sin(100•π•t))/π. Видно, что это незатухающие синусоидальные колебания (благодаря бесконечной добротности).
/*Найдём результат обратного преобразования Лапласа от воздействия ступеньки на резонансный регулятор*/ W1:ilt(V1/p,p,t),Num; /*Подставляем значение для перехода от дискретной передаточной функции к цифровому фильтру*/ V7:factor(V2),z=1/zm; /*Раскрыть скобки у полиномов числителя и знаменателя*/ V8n:expand(num(V7));V8d:expand(denom(V7)); /*Выделить для числителя и знаменателя коэффициенты при степенях zm*/ V9cn:makelist(coeff(V8n,zm,i),i,0,hipow(V8n,zm)); V9cd:makelist(coeff(V8d,zm,i),i,0,hipow(V8d,zm)); /*Разделить числитель и знаменатель на свободное слагаемое знаменателя (при нулевой степени zm) */ V10cn:makelist(i/V9cd[1],i,V9cn); V10cd:makelist(i/V9cd[1],i,V9cd); /*Подставить численные параметры*/ V11N:float(V10cn),Num; V11D:float(V10cd),Num;
Цифровому фильтру будет уделено внимание в отдельной публикации, поэтому кратко можно описать следующий процесс, как показано на рисунке 12.
Переход:
подставить zm = z⁻¹ в дискретную передаточную функцию и упростить
выделить полиномы числителя и знаменателя
сформировать массив коэффициентов при zm начиная с нулевой степени (свободное слагаемое)
разделить коэффициенты числителя и знаменателя на свободное слагаемое знаменателя, тем самым, произведя нормировку цифрового фильтра
Далее следует скрипт, который позволяет реализовать цифровой фильтр в Maxima. Более подробно в следующей публикации.
/*Коэффициенты числителя и знаменателя*/ CnumN:V11N;CdenN:V11D; DigFiltReset(Filter):=block ( [Nx,Ny,i], Cx:Filter[3],Cy:Filter[4], Nx:length(Cx),Ny:length(Cy), for i:1 thru Nx do (Filter[1][i]:0), for i:1 thru Ny do (Filter[2][i]:0), Filter ); DigFiltInit(Cx,Cy):=block ( [Fliter,X,Y,Nx,Ny], Nx:length(Cx),Ny:length(Cy),Filter:makelist(0,4), X:makelist(0,Nx),Y:makelist(0,Ny),Cy[1]:0, Filter[1]:X,Filter[2]:Y,Filter[3]:Cx,Filter[4]:Cy, Filter ); DigFiltStep(Filter,InputX):=block ( [Acc,i,T1,T2,Cx,Cy,X,Y], X:Filter[1],Y:Filter[2], Cx:Filter[3],Cy:Filter[4], Nx:length(Cx),Ny:length(Cy), X[1]:InputX, Cy[1]:0,Acc:0, ( for i:1 thru Nx do (Acc:float(Acc+X[i]*Cx[i])), for i:1 thru Ny do (Acc:float(Acc-Y[i]*Cy[i])), Y[1]:Acc, T1:X[1], for i:1 thru Nx do (T2:X[i],X[i]:T1,T1:T2), T1:Y[1], for i:1 thru Ny do (T2:Y[i],Y[i]:T1,T1:T2) ), Filter[1]:X,Filter[2]:Y, Filter[3]:Cx,Filter[4]:Cy, Filter ); DigFiltCycle(Filter,Function,Niter,Fd):=block ( [Fst,j,Out,OutIn,Fv],Fst:Filter, DigFiltReset(Fst), Out:makelist([0,0],i,Niter),OutIn:makelist([0,0],i,Niter), for j:1 thru Niter do ( Fv:float(ev(subst(t=(j-1)/Fd,Function))), Vhod:1, DigFiltStep(Fst,Vhod), Out[j][2]:Fst[2][1], Out[j][1]:j, OutIn[j][2]:Fv, OutIn[j][1]:j ) , [Out,OutIn] ); /*Инициализация переменных состояния фильтра и коэффициентов*/ Fa:DigFiltInit(CnumN,CdenN); Out:DigFiltCycle(Fa,W1,50,subst(Num,fd))$ wxplot2d([[discrete, Out[1]],[discrete, Out[2]]],[style,points,lines],[y,-10,10],[grid2d,true],[ylabel,"Fp(t),Fz(t)"],[color,green,black],[legend,false]),wxplot_size=[480,320];
В итоге получается требуемый отклик во временной области:
Коэффициент цифрового фильтра 2 cos((2•π•f₀)/fd) можно рассчитывать непосредственно как константу или используя f₀ от блока ФАПЧ.
