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

У

правление подвижными объектами

УДК 517.977.5:629.7

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

А.В. Силаев

Институт проблем управления им. В.А.Трапезникова РАН, г. Москва

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

ВВЕДЕНИЕ

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

Для исследования динамики движения БКК применяется цифровое моделирование. При этом возникает проблема одновременного воспроизве-

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

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

уравнений численными методами происходит накопление погрешности вычислений, что может при выборе слишком большого шага интегрирования привести к мнимой неустойчивости моды или ее огибающей: на самом деле мода (ее огибающая) устойчива, а по результатам моделирования она представляется неустойчивой из-за сильного роста погрешности вычислений. С другой стороны, выбор слишком малого шага приведет к существенному возрастанию времени моделирования движения БКК. Поэтому следует так выбирать шаг интегрирования, чтобы погрешность моделирования на заданном отрезке времени не превышала бы заданного максимально допустимого значения. При этом достигается компромисс между временем моделирования и точностью. Необходимо разработать методику выбора шага интегрирования уравнений, описывающих как сами координаты мод, так и их огибающие, при заданной длительности отрезка времени, на котором проводится моделирование. Данная работа посвящена решению этой задачи.

1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим в качестве математической модели углового движения БКК модально-физическую модель [1], в которой полное угловое движение объекта представляется суммой движения БКК как жесткого тела и упругих мод. В предложенной в работе [2] приближенной модели среди упругих мод выделяются моды, оказывающие существенное влияние на угловое движение объекта в целом, они описываются уравнениями для мгновенных значений их координат. Уравнение движения БКК как жесткого тела вместе с уравнениями существенных мод формируют ядро приближенной модели. Остальные (как правило, высокочастотные) моды, оказывающие менее существенное влияние на угловое движение, описываются своими огибающими и формируют остаточную часть модели.

Таким образом, приближенную модель углового движения БКК можно представить в следующем виде. Ядро модели описывается уравнениями

X = и, (1)

X + ю2 X = к и, / = 1, т, (2)

где х — координата движения БКК как абсолютно жесткого тела, X — координата /-й моды упругих колебаний, и — приведенный управляющий момент, ю; и к{ — соответственно собственная циклическая частота и коэффициент возбудимости /-й моды, т — число мод, входящих в ядро. Предполагается, что собственные частоты мод ю; упоря-

дочены по возрастанию. Отметим, что случай не-

демпфированных мод, рассматриваемый в данной работе, наиболее опасен в смысле влияния на устойчивость всей системы управления ориентацией БКК в целом.

В работе [2] показано, что критерием существенности моды при принятии решения о включении ее в ядро или в остаточную часть может служить так называемая относительная степень возбудимости

моды, равная к./ юг- . Результаты, полученные на основе этого критерия, хорошо согласуются с результатами, достигнутыми при использовании другого общепринятого критерия, основанного на теории внутреннего балансирования [3].

Пусть общее число учитываемых мод равно n (n > m). Тогда оставшиеся n — m мод, менее существенных в смысле их влияния на полное движение объекта, формируют остаточную часть модели и описываются своими огибающими. Возможны различные методы получения этих огибающих. В работе [2] предложены два метода получения огибающей моды упругих колебаний: тригонометрический (при релейно-логическом управлении) и динамического оценивания (при более общем виде закона управления).

Рассмотрим метод динамического оценивания, в соответствии с которым огибающая а. i-й моды определяется следующим образом:

A = —(k;/ra;.)usmra;i, Д- = —(к/га)и cosra.i, (3)

Г~2 ~2 -------

а. = л/A i + Bi , i = m + 1, n, (4)

при начальных условиях A (0) = x (0), Д- (0) = = — X (0)/юг- где начальные условия мод X (0) и X (0) предполагаются заданными.

Движение X (t) ограничено областью [—а .(f), a(t)]. Движение БКК с учетом упругих колебаний описывается уравнением

x = x + X, (5)

где x — измеряемая выходная координата объекта, X — дополнительное изменение координаты, вызванное упругими колебаниями и определяемое в данной приближенной модели суммой координат мод ядра и суммой огибающих мод остаточной части, взятых как со знаком «плюс», так и со знаком «минус»

m n

x = X x ± X a. ()

i = 1 im + 1

Управляющее воздействие u = u(x, X, t) примем релейно-логическим, что обусловлено прежде всего наибольшим распространением релейных исполнительных органов управления БКК, а также

наибольшей опасностью релейного управления в смысле его влияния на амплитуду колебаний. Пусть управляющее воздействие имеет вид следующей функции времени:

N

u(t) = Z Auj1(t — j),

(7)

где 1(t — ]) — единичное ступенчатое воздействие,

подаваемое в момент времени t = ], j = 1, N — номер такта (интервала постоянства) управления. Значения Auj и j до текущего момента времени t (t > tN) заданы.

Необходимо выполнить цифровое моделирование приближенной модели БKK (1)—(7) на заданном отрезке времени [О, T] при условии обеспечения заданной допустимой точности вычислений. Поскольку наибольшая погрешность вычислений накапливается при интегрировании уравнений (2) и (3), то задаются максимально допустимые среднеквадратические отклонения (CKO) вычисленных значений координат мод omax[ X, ], i = 1, m и

огибающих cmax[a;.], i = m + 1, n от их точных значений. Оценки сходимости и погрешности процесса приближенного интегрирования дифференциальных уравнений обычно получаются известными приемами вычислительной математики. Использование CKO погрешности, несмотря на детерминированный характер уравнений, позволяет получить сглаженную оценку погрешности в зависимости от числа шагов интегрирования и размера шага, так как исключает динамику погрешности на малых отрезках времени при сохранении тенденции ее изменения на большом отрезке времени.

Требуется так выбрать шаги интегрирования уравнений (2) и (3), чтобы CKO координат мод и огибающих, найденных с помощью выбранного численного метода, не превышали заданных значений.

2. ВЫБОР ШАГА ИНТЕГРИРОВАНИЯ УРАВНЕНИЯ МОДЫ

При управлении вида (7) аналитическое решение уравнения (2) для /-й моды имеет вид

х-( 0)

х; (^ = х; (0)cosю.t + —— smю. t +

Ю;

ki

N

+ “2 Z Au/1 - cos«i(t - tj)),

юг j = 1

(8)

где tN — последний перед текущим моментом времени t момент переключения управления ^ > tN).

При произвольном виде управления и(^ в выражении (8) последнее слагаемое следует заменить

на

k

ю

■ J smra;(t — x)u(x)dx.

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

Для демонстрации работоспособности предложенной приближенной модели БКК выберем простейший метод численного интегрирования дифференциальных уравнений — метод Эйлера.

Получим приближенное решение уравнения (2) путем численного интегрирования по методу Эйлера. Для этого запишем уравнение (2) в векторной форме:

У (О = 4К0 + Яи(0, (9)

~ ~ Т

где у(^ = (X (^, х (^) — вектор состояния, А =

О 1

2

B =

для простоты не

Здесь и далее индекс /, со-

4 -ю" 0

ответствующий номеру моды, указывается.

Обозначим через к шаг численного интегрирования уравнения моды (2) (или (9)). Тогда в соответствии с методом Эйлера получим

у(к) = (I + Ак)у(0) + Ви(0)к, у(2к) = (I + Ак)2у(0) + (I + Ак)Ви(0)к + Ви(к)к,

y(nh) = (I + Ahfy^) +

n - 1

+ Z (I + Ah)n - 1 - 1Bu(ih)h,

i = 0

(1О)

где I — единичная матрица размером 2x2.

Для простоты рассмотрим сначала случай, когда и(^ = 0. Тогда у(ик) = (I + Ак)”у(0).

Для вычисления у(ик) воспользуемся способом, описанным в книге [4]. Так как матрица А — простая (ее собственные значения 112 = ±/ю), то ее

Т

собственные векторы е12 = (1, ±/ю) образуют базис в пространстве состояний. Поэтому любое начальное состояние у(0) может быть представлено в виде у(0) = а1е1 + а2е2, где а12 = у1(0)/2 + /у2(0)/(2ю). Тогда

2

у(пк) = £ (I + Ак)”а;е;.

o

1

1

Поскольку для собственных векторов е1 и е2 матрицы A справедливо Ae. = 1.е;, то

y(nh) = Z (1 + h) a.e. =

i = 1

= (1 + /кю)”а1 1 1 + (1 — /кю)”а2 1 1.

V / Ю ) V —/ Ю У

После подстановки выражений а1 и а2 для первой координаты вектора состояния получим:

[”/2] • т

уДик) = х(ик) = у1(0) X (-1); С2; (кю)2; +

i = О

+

У2( О)

[(n - 1)/2]

ю

z (—1).с:+1 ^ю)2і + l,

i = О

где [•] — целая часть числа, С^ — число сочетаний

из п элементов по і.

Теперь перейдем к случаю и(/) ^ 0. Каждое слагаемое правой части выражения (10) вида (I + АА/В после аналогичных преобразований можно представить вектором с координатами

к [(1 -1)/2]

Ю

- 1 )/2]

Z (—1). C2 +1 ^ю)2і+1

i = О

[V2] . ,2.

и k Z (—1). C2. ^ю)2г.

і = О

Приближенное значение первой координаты вектора состояния

[ n/2] • 2; т

yl(nh) = Уі(О) Z (—1). сП- ^ю)2г +

i = О

+

У2 (0)

[(n - 1)/2]

ю

Z (—1). C;+1 ^ю)2і + 1 +

i = О

n - 2 [(n - r- 2)/2]

+kkh X u(rh) X (-1)г C«2-+r-1 (hra)2i + 1 .(11)

r = 0 i = 0

Из выражений (8) и (11) получим выражение для погрешности вычисления координаты моды s1(nh) = y1(nh) — x(nh):

,[n/2] -

s1(nh) = x(0) [ X (-1)iCf (hra)2i - cosranhj + i=0

/[(n- 1)/2] N

+ x(0) ^ ^ (-1)- d1 +1 (hra)2- +1 - sinrahn) +

i = О

1 1 n - 2 [(n - r - 2)/2]

+ kh Z u(rh) Z (—1).Cn2-+r-1 ^ю)2і + 1

ю r = О i = О

kN

— —2 Z Auj(1 — cosю(nh — tj)). ю j = l

Для вторых координат векторов у(пА) и є(пА) получены следующие соотношения:

[(«-1)/2] _

У2(пА) = -юуі(0) X (-1)' СП' + 1 (М2г + 1 +

i = О [ n/2]

+

У2(0) Z (—1).сП- ^ю)2і

+

i = О

n - 1 [(n - r - 1)/2]

+ kh Z u(rh) Z (—1). Cn2- r- 1 ^ю)2і, (12)

r = О i = О

,[(n - 1)/2]

/ [("- *-)/ ^] , . 1

£2(nh) = —юХ(0) [ Z (—1). сП- + 1 ^ю)2- + 1 —

i = О

. Лп/2] ... \

- sinюnhJ + X (0) ^ Z (—1). сП- ^ю)2г — cosюnhJ +

i = О

n - 1 [(n - r - 1)/2]

+ kh Z u(rh) Z (—1).Cn2-r-1 ^ю)2і

r = О

i = О

k

ю

N

Z Aujsi^^n — tj).

і = 1

Для объекта, описанного в статье [2], проведено моделирование с целью сравнения координат мод, полученных с помощью выражений (11) и (12), с результатами непосредственного численного интегрирования уравнений (2) по методу Эйлера. Сравнение результатов показало идентичность полученных значений у1(пА), у2(пА) и результатов численного интегрирования.

Для выбора шага интегрирования уравнения моды (2), при котором СКО <г1 = ст[у1] приближенного значения координаты у1 моды на заданном отрезке времени [0, 7] не превышает своего максимально допустимого значения <гшах[ X ], проведем исследование зависимости ст1 от размера шага А и числа шагов п. Исследования проводились на примере объекта, описанного в статье [2], при релейно-логическом законе управления, приведенном там же. В результате было установлено, что при фиксированном значении шага А зависимость ст1 от п близка к линейной до определенного значения п (рис. 1), причем с уменьшением А ширина линейной зоны увеличивается. Это справедливо для различных мод при отсутствии резонанса. При фиксированном значении числа шагов п (рис. 2) значение ст1 с достаточной степенью точности сна-

2 3

чала пропорционально А , затем А и более высоким степеням А. Таким образом, до некоторых критических взаимосвязанных значений п и А, при превышении которых процесс приближенного интегрирования уравнения (2) резко расходится, справедлива следующая аппроксимация: ст1(п, А) « впА ,

2

Рис. 1. Зависимость СКО координаты первой моды от числа шагов интегрирования; h = 0,01 с

Рис. 2. Зависимость СКО координаты первой моды от значения шага интегрирования; п = 2^105

где значение коэффициента р для каждой моды свое. Такая же зависимость наблюдается и для СКО а2 = а[у2] приближенного значения координаты у2 моды, соответствующей скорости X.

Следовательно, при увеличении числа шагов в й раз (й > 0) и одновременном уменьшении шага

в л/й раз значение а1 изменится незначительно. На

этом основана предлагаемая методика выбора шага интегрирования моды:

— по имеющейся зависимости а1(к) при фиксированном числе шагов и* выбирается значение шага к* такое, что а1(к*) < атах;

— вычисляются значения Т* = и*к* и й = = (Т/Т *)3;

— вычисляются число шагов и = и* и

размер шага к = к*/ Цй.

Нетрудно видеть, что при таком выборе и и к и справедливости зависимости а1(и, к) « рик получим а1(и, к) « а1(и*, к*) < атах и ик = Т.

3. ИССЛЕДОВАНИЕ ПОГРЕШНОСТИ ВЫЧИСЛЕНИЯ ОГИБАЮЩЕЙ МОДЫ

В соответствии с методом динамического оценивания огибающая моды получается путем интегрирования уравнений (3) при управлении (7). При этом аналитическое решение имеет вид

к Г N — 1

А (0 = А (0) + --2 V X и^еоБю^. + 1 - еоБю^.) +

ю У = 1

+ uN(cosюt — СОБЮ^)), (13)

к ( N — 1

В (0 = В (0) — — ^ X и^тю. + 1 — sinюtj.) +

ю У = 1

+ uN(sinюt — sinюtN) ), (14)

где / — номер такта управления, / = 1, N; . — момент переключения управления, = 0, t > ^; и. — значение управления на /-м такте.

Приближенное решение системы (3), полученное по методу Эйлера, в момент времени t = ик (к — шаг интегрирования) имеет вид:

кк п — 1

^(ик) = £.(0) — — X и(/к^тю/к,

1 1 ю ^

kh

i = О n - 1

z2(nh) = z2(0) — — Z ^і^шзюй,

ю

і = 0

где г1 и £2 — приближенные значения А и В соответственно; г1(0) = А (0), г2(0) = ВВ (0).

Отсюда получим выражения для погрешностей с а (пА) = г1(пА) — А (пА) и с в (пА) = г2(пА) — В (пА)

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

кА п -1

с а (пА) = — — X и(/А^тю/А —

Ю і = 0

к г 1

""2 І X И/(^ю$ + 1 — ^ю$) +

ю

j = 1

+ uN(cosюt — cosюtN) J,

(15)

кк п 1

с ь (ик) = — — X u(/к)cosю/к +

- ю ^

; = о

к ( N—1

+ — I X и^Шю . + 1 — sinюtj.) + ю У = 1

+ uN(sinюt — sinюtN)), (16)

где ик > tN.

Погрешность са вычисления огибающей я(0 моды, определяемой выражением (4), может быть вычислена следующим образом:

Ев(ик) - Н 8 - (ик) + ^ 8 - (ик) =

а 5А - 5,8 -

А с-( ик) + ,8с -(ик)

л/А2 + В2

(17)

где значения А и В вычисляются в момент времени t = ик. Таким образом, соотношения (13)—(17) позволяют оценить указанную погрешность в момент времени, кратный шагу интегрирования к.

Для выбора шага интегрирования на заданном отрезке времени, при котором обеспечивалось бы допустимое значение СКО огибающей, необходимо исследовать влияние размера шага к и числа шагов и на СКО огибающей а[а(ик)].

Для различных мод объекта [2] получены следующие результаты. При фиксированном шаге к зависимость а(и) носит нелинейный и немонотонный характер с тенденцией к возрастанию (см., например, зависимость а(и) для четвертой моды при к = 0,1 с, представленную на рис. 3). При фиксированном числе шагов и зависимость а(к) для всех мод объекта при отсутствии резонанса носит качественно схожий характер (рис. 4): величина а возрастает медленно, при некоторых значениях шага к, равных для первой моды 0,495; 0,560 с и ряду других значений, имеет место резкий пик а. Это позволяет дать рекомендации для выбора значения шага интегрирования и числа шагов.

4. ОЦЕНКА СНИЖЕНИЯ ВРЕМЕНИ МОДЕЛИРОВАНИЯ

Оценим преимущество предложенной приближенной модели, выражающееся в снижении времени моделирования движения БКК на заданном отрезке времени. Пусть т0 и т — соответственно время моделирования одной моды и ее огибающей путем численного интегрирования уравнений (2) и (3) на заданном отрезке времени [0, Т ]. Предположим, что моделирование всех т мод ядра модели осуществляется с одинаковым шагом к0, а моделирование огибающих всех и — т мод остаточной части — с шагом к. Тогда время моделирования

Рис. 3. Зависмость СКО огибающей четвертой моды от числа шагов интегрирования; h = 0,1 с

Рис. 4. Зависимость СКО огибающей первой моды от значения шага интегрирования; п = 105

приближенной модели равно тт0 + (и — т)т. Сравним его со временем моделирования точной модели, полученной путем численного интегрирования уравнений (2) для всех и мод с шагом к0, которое равно ит0.

Будем количественно оценивать преимущество предложенной модели отношением п указанных времен:

_ т т 0 + ( и — т) т = т + (1 — т^

П = > - , ■

и Т 0 и 1 иУ Т 0

Очевидно, что отношение т/т0 обратно пропорционально отношению шагов интегрирования

уравнений огибающей (3) и моды (2), т. е. — = а ,

т0 к

где а показывает, во сколько раз время моделирования огибающей моды на одном шаге интегрирования больше времени моделирования самой моды на одном шаге. В результате исследований было установлено, что при получении огибающей моды по методу динамического оценивания (т. е. по уравнениям (3) и (4)) а - 5 [2]. Таким образом,

_ т , (л т^| к0

П - - + I 1 — —)а-т .

и 1 иУ к

Отсюда нетрудно видеть, что приближенная модель обеспечивает снижение времени моделирования упругой составляющей движения объекта (т. е. п < 1) при к > ак0.

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

ЗАКЛЮЧЕНИЕ

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

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

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

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

Экспериментально установлена приблизительная зависимость СКО вычисленного значения координаты моды от числа шагов и размера шага интегрирования, позволяющая путем исследования данного СКО на относительно малом отрезке времени при некотором значении шага интегрирования так выбрать шаг, чтобы СКО координаты моды на заданном (большом) отрезке времени не превышало бы максимально допустимого значения.

ЛИТЕРАТУРА

1. Рутковский В.Ю., Суханов В.М. Большие космические конструкции: модели, методы исследования и принципы управления // Автоматика и телемеханика. — 1996. — № 8. — С. 55—66.

2. Силаев А.В., Суханов В.М. Методы оценки огибающей упругих колебаний деформируемого космического аппарата // Там же. — 2002. — № 4. — С. 76—90.

3. Грегори-мл. К.З. Использование теории внутреннего балансирования для редукции моделей больших упругих космических летательных аппаратов // Аэрокосмическая техника. — 1985. — Т. 3. — № 6. — С. 71—80.

4. Бахвалов Н.С. Численные методы. — М.: Наука, 1975.

® (495) 334-87-79, e-mail: [email protected]

Статья представлена к публикации членом редколлегии

В.Ю. Рутковским. □

V Не заІЇцуьіїе подписаться!

1 л V у ч Подписку на журнал "Проблемы управления" можно оформить с любого месяца в любом почтовом отде-

1 | лении (подписные индексы 80508 и 81708 в каталоге Роспечати или 38006 в объединенном каталоге "Пресса Рос-

ю. и н а в о б е р е ж у м о в р е о е а ы с ы в я и ц к а O е р а р е м о н е ы н ь ел O ит ю. и ц к а O е р ез ер ч е ж к а а "), и" и с

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

можно приобрести любой номер журнала за прошедшие годы. Позвоните в редакцию по тел. (495) 334-92-00 или пришлите Ваш

адрес по электронной почте [email protected] - - и подписка будет оформлена за один рабочий день. Почтовые расходы по пересылке

журнала редакция берет на себя.

Оформить подписку можно и через распространяющие наш журнал во многих городах России и СНГ агентства альтер-

нативной подписки:

Альтерра-информ (г. Казань) Коммерсант-Курьер "ПрессМоушн"

в (843) 292-87-09, 512-24-45 www.komcur.ru www.pressmotion.ru

e-mail:[email protected] Головной офис в г. Казань в (495) 221-08-77 (доб. 100)

Артос-гал в (8432) 91-09-82 e-mail: [email protected]

www.setbook.ru Красносельское агентство Урал-Пресс

в (495) 995-44-63, 504-13-45 "Союзпечать" (г. Екатеринбург)

Вся пресса (Корпоративная подписка) www.ural-press.ru

в (495) 234-03-08 в (495) 783-29-60 в (343) 262-61-35,

Интер-почта Орикон-М 262-65-43

www.interpochta.ru в (495) 937-49-59, 937-49-58 Экспресс

в (495) 500-00-60 e-mail:[email protected] в (495) 234-23-80,

Пресс Инфо Орикон-С (г. Саратов) 783-90-29

в (495) 256-26-06 в (8452) 37-67-87, (8452) 51-83-65 e-mail: [email protected]

e-mail: [email protected] 1