-
Notifications
You must be signed in to change notification settings - Fork 0
/
thesis.tex
executable file
·4365 lines (3560 loc) · 341 KB
/
thesis.tex
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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
\documentclass[12pt,lot, lof]{puthesis}
\title{Modeling of plasma rotation control for NSTX and NSTX-U}
\submitted{September 2016} % degree conferral date (January, April, June, September, or November)
\copyrightyear{2016} % year in which the copyright is secured by publication of the dissertation.
\author{Im\`ene Rym Goumiri}
\adviser{Clarence W. Rowley and David A. Gates} %replace with the full name of your adviser
%\departmentprefix{Program in} % defaults to "Department of", but programs need to change this.
\department{Mechanical and Aerospace Engineering}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\
%%%% Tweak float placements
% From: http://mintaka.sdsu.edu/GF/bibliog/latex/floats.html "Controlling LaTeX Floats"
% and based on: http://www.tex.ac.uk/cgi-bin/texfaq2html?label=floats
% LaTeX defaults listed at: http://people.cs.uu.nl/piet/floats/node1.html
% Alter some LaTeX defaults for better treatment of figures:
% See p.105 of "TeX Unbound" for suggested values.
% See pp. 199-200 of Lamport's "LaTeX" book for details.
% General parameters, for ALL pages:
\newcommand{\vect}[1]{\mathbf{#1}}
\newcommand{\R}{\mathbb{R}}
\renewcommand{\topfraction}{0.85} % max fraction of floats at top
\renewcommand{\bottomfraction}{0.6} % max fraction of floats at bottom
% Parameters for TEXT pages (not float pages):
\setcounter{topnumber}{2}
\setcounter{bottomnumber}{2}
\setcounter{totalnumber}{4} % 2 may work better
\setcounter{dbltopnumber}{2} % for 2-column pages
\renewcommand{\dbltopfraction}{0.66} % fit big float above 2-col. text
\renewcommand{\textfraction}{0.15} % allow minimal text w. figs
% Parameters for FLOAT pages (not text pages):
\renewcommand{\floatpagefraction}{0.66} % require fuller float pages
% N.B.: floatpagefraction MUST be less than topfraction !!
\renewcommand{\dblfloatpagefraction}{0.66} % require fuller float pages
% The documentclass already sets parameters to make a high penalty for widows and orphans.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\
%%%% Use packages
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage{tikz}
\usepackage{adjustbox}
\usepackage{overpic}
\usepackage[numbers]{natbib}
%\usepackage{booktabs}
\graphicspath{{figs/}{figs/chap7/}}
\usetikzlibrary{shapes,arrows,positioning,fit,calc}
\tikzstyle{block} = [draw, rectangle, fill=yellow!20,
minimum height=2em, minimum width=3em, inner sep=0.5em]
\tikzstyle{sum} = [draw, circle, inner sep=0pt, minimum size=0.6em]
\tikzstyle{junction} = [draw, circle, fill, inner sep=0pt, minimum size=1mm]
% Uncomment the following to leave junctions "plain" (no small filled in circle):
\tikzstyle{coord} = [coordinate]
\tikzstyle{connector} = [->,thick]
\tikzstyle{tag} = [node distance=1mm]
\newcommand{\remark}[1]{\textcolor{blue}{[#1]}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Printed vs. online formatting
\ifdefined\printmode
% Printed copy
% url package understands urls (with proper line-breaks) without hyperlinking them
\usepackage{url}
\else
\usepackage{hyperref}
\hypersetup{bookmarksnumbered,colorlinks,allcolors={black},urlcolor={blue}}
\makeatletter
\hypersetup{pdftitle=\@title,pdfauthor=\@author}
\makeatother
\fi % printed or online formatting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\
%%%% Define commands
% Define any custom commands that you want to use.
% For example, highlight notes for future edits to the thesis
%\newcommand{\todo}[1]{\textbf{\emph{TODO:}#1}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\
%%%% Front-matter
% For early drafts, you may want to disable some of the frontmatter. Simply change this to "\ifodd 1" to do so.
\ifodd 0
% front-matter disabled while writing chapters
\renewcommand{\maketitlepage}{}
\renewcommand*{\makecopyrightpage}{}
\renewcommand*{\makeabstract}{}
% you can just skip the \acknowledgements and \dedication commands to leave out these sections.
\else
\abstract{
This thesis studies some applications of feedback control design for plasma physics, specifically plasma drift waves and plasma toroidal rotation and covers two major steps: reduced-order modeling and controller design.
This dissertation focuses mainly on the toroidal plasma rotation but begins with the Hasegawa-Wakatani (HW) problem, a classic model of plasma drift waves zonal flows coupling, as a preliminary case study where the basic methodology of reduced-order modeling and control is applied to demonstrate the effectiveness and applicability of the approach.
First, the development of a model-based feedback control that stabilizes an unstable equilibrium in the HW equations is studied: a balanced truncation (a model reduction technique) is applied to obtain a low-dimensional model of the linearized HW equations. Then a model-based feedback controller is designed for the reduced order model using a Linear Quadratic Estimator (LQE) which only requires a small set of sensors. Results show that this controller applied to the original non-reduced nonlinear HW equations stabilizes the equilibrium and suppresses the transition to drift-waves instabilities.
Then the thesis dives into the core subject which is the control of plasma toroidal rotation in tokamaks. It uses experimental measurements from the National Spherical Torus Experiment (NSTX) and is aimed at controlling plasma rotation using two different types of actuation: momentum from injected neutral beams and neoclassical toroidal viscosity generated by three-dimensional applied magnetic fields. Based on the data-driven model obtained, a feedback controller is designed, and predictive simulations using the TRANSP plasma transport code show that the controller is able to attain desired plasma rotation profiles given practical constraints on the actuators and the available measurements of rotation.
The last part studies the rotation control on the upgraded device NSTX-U. The major change comes from the addition of a second neutral beam injector which adds three more actuators to the designed controller and thus gives us considerably more flexibility, at the expense of added complexity in the modeling and control of simultaneously the toroidal rotation and the stored energy.
Because NSTX-U modeling is a model-based design (we rely heavily on model predictions and sensors measurements) and experimental data from NSTX-U are not available, a study of the robustness of our controller to some parameters uncertainties in particular the perpendicular momentum diffusivity profile $\chi_{\phi}$ and the confinement time $\tau_e$ is developed.
}
\acknowledgements{
First and foremost, I would like to thank my advisors Professors Clarence W. Rowley and David A. Gates for their constant guidance and support over these past six years. Whenever I needed help or had a question, they were always here for me, spending time discussing with me research problems in details. I am very grateful for the constant encouragement Clancy and David have provided; it helped me believe in myself, and motivated me to push my limits further. I am also very grateful for their patient teaching on coding, how to write papers, collaboration, presentations and talks, and even Matlab and \LaTeX\ skills.\\
%
I would like to thank my readers Prof.\ Egemen Kolemen and Dr.\ Jonathan E. Menard for their helpful comments and Prof.\ Naomi E. Leonard and Dr.\ Stefan P. Gerhardt for being my examiners.\\
%
I would like to thank Prof.\ Steve A. Sabbagh for his valuable and constant help, comments and suggestions throughout my thesis projects. His support and involvement made him a very important mentor. I want to thank Prof. N. Jeremy Kasdin and Dr.\ Stefan P. Gerhardt for being involved in my committee and following my progress closely during these years. I am also thankful to Prof.\ John A. Krommes and Dr.\ Jeffrey B. Parker for their help and time on the Hasegawa-Wakatani project. A special thanks goes to Dr.\ Mark D. Boyer for his precious advices, feedbacks and readings; our collaboration was very instructive and I hope will carry on beyond this work.\\
%
I would like to thank the amazing members of Prof. Rowley's lab. In particular, I would like to thank Dr.\ Brandt Belson and Dr.\ Lauren E. Padilla not only for many helpful and inspirational discussions but also for sharing up and down moments over theses years. I would like to thank Dr.\ Zhanhua Ma, Prof.\ Steven Brunton, Dr.\ Matthew Williams, Prof.\ Mark Luchtenburg, Prof.\ Maziar Hemati, Dr.\ Kevin Chen, Dr.\ Jonathan Tu, Scott Dawson for their friendship and good mood. I would like to thank Anthony DeGennaro, a dear friend who helped me with advices and talking and even some driving during tough times. \\
%
Special thanks for my graduate administrator Jill F. Ray for her tremendous help and support at anytime of the day. She really brings joy and positiveness to the MAE department.\\
%
I also want to thank my family for their loving support, especially to my parents for their encouragement from Algeria. To my mom and dad, I am here today thanks to your unflagging emphasis on education over almost three decades. Finally, with all my heart, I would like to thank my husband Simon, who has steadfastly walked every step of this endeavor with me; thank you for your constant support, friendship, and love. I would never be where I am in my life without you.\\
%
A final word of thanks goes to all of those who helped me in any capacity while I completed this work, especially Dr. Robert Andre and Dr. Robert Budney who were my reference whenever I had to run my code in TRANSP. Thank you for your help with coding and debugging.\\\\
%
This dissertation carries the number 3311T in the records of the Department of Mechanical and Aerospace Engineering.
}
\dedication{To Simon, Lila and Malik. \\ To my parents.}
\fi % disable frontmatter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\
%%%% Hide some chapters
%%% If you want to produce a pdf that includes only certain chapters, specify them with includeonly, in addition to including all chapters below.
%\includeonly{ch-intro/chapter-intro}
%%% You can also specify multiple chapters.
%\includeonly{ch-intro/chapter-intro,ch-usage/chapter-usage}
%\includeonly{chap1,chap2,chap3}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Notes:
% Footnotes should be placed after punctuation.\footnote{place here.}
% Generally, place citations before the period~\cite{anotherauthor}.
% The proper usage for i.e., and e.g., include commas ``(e.g., option A, option B)''
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Import chapters
\begin{document}
\makefrontmatter
% If you've disabled frontmatter, you can insert the toc manually
%\tableofcontents\clearpage
% \include lets us split up the document (and each include starts a new page):
%\include{ch-intro/chapter-intro}
\part{Plasma modeling and control}
\chapter{Introduction}
\section{Fusion and Tokamaks}
It became general knowledge that the fossil fuel era will end during this century. Coal, petroleum and natural gas are becoming harder and deeper to extract. Scientists agree on the fact that a shortage of fossil fuel energy is inevitable in less than 30 years from now as shown in figure \ref{projected_energy_shortfall}~\cite{GIEC}.
\begin{figure}
\centering
\includegraphics[width=\linewidth]{projected_energy_shortfall.pdf}
\caption{Planned energy shortfall. (Figure courtesy of Lawrence Livermore National Laboratory. \cite{GIEC})}
\label{projected_energy_shortfall}
\end{figure}
Renewable energy sources, such as solar and wind power are well ranked to contribute to energy needs in both mature and emerging technologies, but energy from these sources is unstable, intermittent and insufficient to satisfy our world increasing greed of power.
Nuclear fission and fusion, two opposite reactions that occur by splitting heavy atoms such as uranium or fusing light ones such as hydrogen respectively, can potentially produce enough energy to supply the world power demands.
Although fission is a mature technology that is commercially available in all the nuclear plants, its by-products are highly radioactive and long lasting and its reactors are at risk of nuclear accidents due to large uncontrolled release of energy (Tchernobyl 1986, Fukushima 2011). Finally, the main nuclear fuel for fission, Uranium-235 or Plutonium-239 are nonrenewable resources (like all metals), non-recyclable, and not always easily exploited under economical and environmental acceptable conditions.
Fusion, on the other hand, offers undeniable advantages, and is the best candidate for an abundant and clean source of energy: no environmental pollution during operations since it releases helium, no risk of a nuclear accident since uncontrolled release of energy seen in fission cannot happen and finally a sufficient nuclear fuel supply due to the use of deuterium and lithium, both extremely abundant and easy to extract from nature. All these arguments enables us to state that fusion is the quickest energy solution available to us.
The most popular device where the fusion process occurs is called \textbf{Tokamak} and this is exactly what we will focus on for the rest of this dissertation.
Nuclei do not naturally fuse. They are positively charged so they repel each other. In order to fuse, Nuclei have to overcome this huge electrostatic repulsion before they can get close enough together such that the strong nuclear force which maintains nuclei together can kick in.
In order to overcome the natural electrostatic barrier, we have to make the fuel nuclei six to seven times hotter than the Sun's core. In fact to initiate any fusion reaction, a temperature of about 120 million degrees Celsius is required. At such extremely high temperatures, the fuel atoms are dissociated into their component electrons and nuclei, forming the fourth state of matter called \textbf{plasma}.
Keeping this plasma in one place long enough for the nuclei to fuse together is not trivial. Therefore the main idea of the tokamak device is to confine plasma using strong magnetic fields, generated by coils of electrical superconductors around a donut-shaped magnetic bottle in which the plasma is trapped: the combined effect of electric current flowing in the toroidal and poloidal field coils and in the plasma itself produces helical magnetic field lines that create a configuration from which plasma charged particles never leave the torus (inside the tokamak).
Plasmas in tokamak devices need to be confined in order to have a better fusion reaction so increasing the confinement time is a key feature as the longer the plasma is confined, the more fusion energy is extracted.
All existing tokamaks are pulsed devices which means that the plasma is maintained within the tokamak for only a few seconds to several hours (TRIAM-1M maintained a plasma discharge of over 100\,kA for over 5 hours using lower hybrid current drive). Each of these operating cycles of heating the plasma that cools down is called a \textbf{discharge} or a \textbf{shot}. One goal is to increase the duration of the shot. Today's plasma experiments in tokamaks can confine plasmas at the required temperatures for net power gain, but the plasma density and energy confinement time (a measure of the cooling time of the plasma) are too low for the plasma to be self-heated (a burning plasma).
In current fusion experiments, neutral beams, which consist of uncharged atoms of deuterium at high velocity, are being used to heat the plasma. These particles collide with the moving particles of the plasma transferring their momentum to the background plasma and eventually further heating the plasma. This is the most frequently used method of heating the plasma \cite{Stix72, Wagner82, Kaye85, Strachan87}. In fact this is an important element that we will be using in our rotation research, not only as a heater for the plasma but an important contributor of plasma rotation as during collisions, neutral beams also transfer momentum \cite{Suckewer79, Suckewer81, Strait07}. Others methods like ohmic heating through induction, or radio-frequency (RF) heating can also be used for plasma heating. We don't expand on these heating methods but more details can be found in \cite{Stix75, Fisch78, Kubo83, Goldston84}.
It is very important for safety and control to be able to monitor the different plasma properties during a shot, so we can adapt our response in real-time, and understand better how the distributions of the plasma properties behave during the experiment.
Plasma temperature is so hot inside the tokamak. Putting any measurement device during an operation can cause its evaporation or melting.
Luckily, we can rely on indirect methods of measurement and diagnostics, and we can also take advantage of our knowledge of plasma axisymmetric force balance which gives rise to magnetic flux contours and where pressure is constant. Controlling the magnetic flux enables us to control the plasma.
The magnetic reconstruction of these inner flux contours is done using external probes measuring the localized magnetic fields and current flows. However these external measurements are not sufficient to deduce all internal parameters so they have to be complemented by measurements made using different methods. The two most frequent methods are using lasers to see how light is scattered or slowed as it passes through the plasma or sending a beam of neutral atoms through the plasma and analyzing the resulting optical emissions.
This allows us to measure temperature, pressure, density, magnetic field and current throughout the plasma. See \cite{Equipe78, Orlinskij88, Matthews94, McKee99, Hutchinson02} for more details.
\section{NSTX and NSTX-U devices}
A spherical tokamak is a type of fusion power device based on the tokamak principle. It is notable for its very low \textbf{aspect ratio}. A traditional tokamak has a toroidal confinement area that gives it an overall shape similar to a donut with a large hole in the middle. The spherical tokamak reduces the size of the hole to almost zero, resulting in a plasma shape that is almost spherical, often compared to a cored apple. The spherical tokamak is sometimes referred to as a spherical torus and often shortened to \emph{ST}. A comparison against a conventional tokamak is shown in Figure~\ref{nstx1}.
\begin{figure}[htbp]
\centering
\begin{tikzpicture}
\node[anchor=south west,inner sep=0] at (0,0) {\includegraphics[width= 0.7\linewidth]{nstx}};
\draw[red, very thick] (5.47,1.7) -- ++(down:2);
\draw[red, very thick] (6.67,1.7) -- ++(down:2);
\draw[red, very thick] (7.52,1.7) -- ++(down:2);
\draw[red, thick, <->] (5.47,-.2) -- node[black, fill=white, inner sep=2pt] {$R$} (6.67,-.2);
\draw[red, thick, <->] (6.67,-.2) -- node[black, fill=white, inner sep=2pt] {$a$} (7.52,-.2);
\end{tikzpicture}
\caption{Sperical tokamak vs conventional tokamak. (Figure courtesy of Culham Centre for Fusion Energy.)}
\label{nstx1}
\end{figure}
The National Spherical Torus Experiment (NSTX) (Fig~\ref{nstx2}) is an innovative magnetic fusion device based on the spherical tokamak concept. It was constructed by the Princeton Plasma Physics Laboratory (PPPL) in collaboration with the Oak Ridge National Laboratory, Columbia University, and the University of Washington at Seattle.
\begin{figure}[htbp]
\centering
\includegraphics[width= 0.9\linewidth]{nstx2}
\caption{Picture of the National Spherical Torus Experiment (NSTX). (Courtesy of PPPL.)}
\label{nstx2}
\end{figure}
NSTX has a low aspect ratio of $A= R/a =1.31$, with the major radius $R =0.85$\,m and the minor radius $a =0.65$\,m.
The experimental NSTX device has several advantages including plasma stability through improved confinement, but requires a very careful design of the toroidal and poloidal field coils, vacuum vessels and plasma-facing components. Moreover, this innovative plasma configuration has the advantage of being able to confine a higher pressure plasma than a conventional tokamak of high aspect ratio for a given confinement magnetic field strength. Since the amount of fusion power produced is proportional to the square of the plasma pressure, the use of spherically shaped plasmas could allow the development of smaller, more economical and more stable fusion reactors. NSTX's attractiveness may be further enhanced by its ability to reach a high \emph{bootstrap} electric current. This self-driven internal plasma current would reduce the power requirements of externally driven plasma currents required to heat and confine the plasma.
The U.S. Department of Energy's Princeton Plasma Physics Laboratory (PPPL) runs the National Spherical Torus Experiment (NSTX), which has undergone a \$94 million upgrade which was completed in 2015 which makes it the most powerful spherical tokamak in the world. Experiments are currently testing the ability of the upgraded spherical facility to maintain a high-performance plasma under extreme heat and power conditions. Results could strongly influence the design of future fusion reactors.
\begin{figure}[htbp]
\centering
\includegraphics[width= 0.85\linewidth]{nstx3}
\caption{NSTX-U cross section (Figure courtesy of PPPL.)}
\label{nstx3}
\end{figure}
The primary components of the upgrade are the complete replacement of the center stack (Figure~\ref{nstx3}), (which consists of 36 22-footlong, 350-pound copper conductors which comprise the inner-leg of the toroidal field (TF) coils, the Ohmic heating (OH) solenoid, and some divertor coils), and the addition of a second neutral beam injector (Fig~\ref{nstxu}), aimed more tangentially compared to the present original set of beams, and this will give us considerably more flexibility and power to heat the tokamak.
%Temperature inside the NSTX-U exceeds the 15 million degree Celsius core of the sun and NSTX-U magnets are 20,000 times stronger than the Earth's magnetic field.
\begin{figure}[htbp]
\centering
\includegraphics[width= 0.6\linewidth]{NSTXU}
\caption{Cut-away top view of NSTX-U showing the trajectory of the neutral beams. (Figure courtesy of PPPL.)}
\label{nstxu}
\end{figure}
\section{Various control problems in Tokamaks}
During a tokamak pulse, several things happen:
\begin{itemize}
\item A plasma is created.
\item The plasma is ramped up to the reference flat top current and heated to the ignition temperature needed.
\item The plasma is cooled down (\emph{ramped down}) and terminated.
\end{itemize}
The key step of maintaining the plasma in a current flat-top is not quite reached yet \cite{Green03}, but all the other steps are mastered. During a pulse, time variation of coil currents and plasma geometrical and physical parameters are the predetermined nominal values computed from magnetohydrodynamic (MHD) simulation codes to determine the magnetic field and plasma current density necessary for the equilibria \cite{Lao85, Takeda91, Lao05}.\\\\
%
The performance objectives in any tokamak are:
\begin{itemize}
\item High plasma temperature
\item High plasma pressure $\beta$
\item Long energy confinement time $\tau_E$
\item High driving plasma current
\end{itemize}
and all these quantities in steady state configuration. Therefore feedback control of the basic functions of plasma initiation, shaping, heating current drive, stabilization and safe termination of discharges becomes necessary especially because tokamaks often operate near the stability limits. Control is also highly needed due to uncertainty and complexity in the mathematical models of the plasma dynamics and unpredictable disturbances, and this is what makes the overall control problems very challenging and very interesting; the most desirable high performance regimes from the perspective of a fusion power reactor tend to be those that are the nearest to instability.
There have been huge accomplishments and advances in reduced-order modeling and control for plasma fusion on different tokamak devices, including vertical position control, shape control, kinetic and current profile control, MHD (magnetohydrodynamic) stabilization and plasma transport reduction.
These control engineering problems can be classified into two groups:
\begin{itemize}
\item \textbf{Electromagnetic control:} which refers to controlling the magnetic and electric fields to regulate the position and shape of the plasma, as well as the total plasma current.
\item \textbf{Kinetic control:} which refers to controlling fueling rates and auxiliary heating to modify the plasma density, temperature, and pressure.
\end{itemize}
This section addresses these different control problems and gives an short overview of the methods used on different tokamaks for each control problem.
\subsection{Electromagnetic control}
\subsubsection{NTM control}
In a high pressure resistive plasma, nested magnetic surfaces can go unstable by tearing and reconnecting with each others. This creates magnetic islands which connect hotter core regions of the plasma to colder regions. This phenomena known as the Neoclassical Tearing Mode NTM is a short circuit effect that drives the plasma heat to leak out of the core, creating a flattening in plasma temperature, pressure and current profile and results in bad confinement and plasma disruption.
%NTM can be caused by seed islands growing from ELMs instabilities or by background plasma turbulence, which can be impossible to eliminate.
The fundamental idea of NTM suppression comes from experimental observations and it is done through active control laws that use highly localized Electron Cyclotron Current Drive (ECCD) to compensate for the plasma current loss. Restoring the current drive would shrink these magnetic islands, restore the nested flux surfaces and stabilize the mode.
In DIII-D for example, a plasma control system based on a ‘‘search and suppress’’ algorithm (active tracking and target lock routines) was able to make either small rigid radial position shifts of the entire plasma or small changes in the toroidal field to find and lock the optimum position for complete island suppression \cite{Haye02}.
%This was based on real-time measurements and the experimental tests showed the first use of active feedback control to provide continuous precise positioning for NTM control.
While NTM control in tokamaks is still in constant progress, it benefited from advances in control algorithms, estimation, real-time computation, actuator technology and diagnostic signal interpretation \cite{Bongers09, Volpe09, Zohm07}.
Another example, in the ASDEX Upgrade device \cite{Rapson13, Rapson14}, the NTM dynamical system (a second order damped system) was simulated in SIMULINK and an anti-windup Proportional-Integral (PI) controller was designed and tuned (optimized) through various simulation in order to obtain the best system performances.
The controller was designed such that it tries to generate invalid input commands, the mode of operation switches to a fallback mode which requires no inputs and is safe by design. It returns to normal functionality if all input commands are valid again. This functionality which is included as part of a standard feature within the ASDEX Upgrade control system bears resemblance with how we handle invalid input commands in our rotation controller where the inputs (beam powers and coil current) cannot exceed a certain physical threshold.
\subsubsection{ELM control}
An Edge Localized Mode (ELM) is a disruptive instability occurring in the edge region of a tokamak plasma when it is heated above a certain power level (high energy confinement, H-mode). ELMs can cause a loss of up to 10\% of the total plasma stored energy on a very short timescale. Therefore, the control of ELMs is an important issue to consider.
A standard ELM control technique consists of imposing a high critical threshold of pressure gradient to the plasma, and keep the input power low so that this critical threshold is not reached.
The drawback of this technique is that ELMs can be beneficial as they naturally reduce the sources of impurities, thus totally suppressing them would provoke an accumulation of impurities that can deteriorate the plasma confinement.
An ideal control strategy would allow to eliminate the ELMs while providing an alternative method for reducing these impurities through density control. There are no sufficiently advanced methods using a feedback control approach yet but many experimental control approaches that have been tested. For example, injecting deuterium pellet for the ASDEX-U or DIII-D tokamaks \cite{Kallenbach05, Loarte14} or controlling by 3D edge magnetic field perturbations \cite{ Loarte14} all allow a significant reduction of ELMs.
\subsubsection{RWM control}
The Resistive Wall Mode (RWM) is one of the major tokamak non-axisymmetric instabilities.
Active feedback of RWMs in tokamaks was investigated in \cite{Fransson00} using control
theory. Control systems were designed to stabilize the resistive wall mode and meet certain performance specifications for a set of test equilibria. A control problem for $n > 0$ RWMs in tokamaks (through system identification) was formulated then several different controllers: P (proportional), PD (proportional-derivative) and $H_{\infty}$ were applied.
\cite{Garofalo01} for DIII-D and \cite{Sabbagh04} for NSTX applied similar active feedback control of RWM approaches based on simulation and experimental tuning of the controller gains to meet some performance requirements.
A combined algorithm for resistive wall mode identification using both a matched filter and a Kalman filter was implemented in the DIII-D plasma control system (PCS) \cite{In06}. The Kalman filter was based on an eigenmode approach which is similar to our model-based reduced-order control approach, it enables to build a low dimensional controller with an observer takes advantage of the relevant magnetic sensors used for measurements.
\subsubsection{Vertical position control}
It has been proven theoretically that plasmas with a vertically elongated cross section improve and increase energy confinement time, and the first experiments performed in the seventies on tokamaks confirmed these theoretical predictions \cite{Pironti05,Robinson78}. A vertical elongated plasma implies a vertical axisymmetric MHD instability that can be stabilized by surrounding the plasma with a superconductive wall \cite{Mori87, Okabayashi74, Blum81} and applying an active feedback system \cite{Rebhan78}. A detailed investigation of techniques for active stabilization of elongated plasmas can be found in Humphreys and Hutchinson's work on the Alcator C-Mod device \cite{Humphreys89}.
The problem of plasma vertical position stabilization and shape control under actuation saturation in the DIII-D Tokamak at General Atomics was studied in \cite{Schuster03} where an anti-windup compensator was designed for a given predesigned nominal plasma vertical position controller guaranteeing global vertical stabilization of the plasma in the presence of actuator saturation for all reference commands.
\subsubsection{Current profile control}
One of the first work on current control has been the work of Gran \cite{Gran77} on the TFTR device (the ancestor of NSTX), where a control system design has been developed using linear optimal control techniques. On that work, both the ohmic heating and equilibrium field coils were controlled to maintain plasma current and plasma position at their desired values.
Due to its effect on confinement, plasma stability, and non-inductively driven plasma current, the control of current profile became critical and a huge advancement in controlling this profile has recently been made on several machines.
For JET \cite{Moreau03, Moreau08, Laborde05}, a system identification procedure has been developed, a system discretization was therefore performed through an expansion onto a finite set of appropriate basis functions and a Galerkin scheme and a state space model structure was obtained. A combining control law of a slow (proportional + integral) and a fast (proportional) feedback loop was then designed to reach the closest self consistent achievable state defined by the minimization of a quadratic integral error.
For DIII-D, \cite{ Boyer13} has designed a current profile controller using a first-principles-driven dynamic model (with minimal parameters determined from experiment). The feedback controller was designed to complement any arbitrary set of feedforward inputs and drive the spatial profile to the desired target profile (reference tracking). Through a nonlinear transformation of the inputs and spatial discretization, a finite dimensional, time-varying linear model for the profile error was obtained. A singular-value decomposition (SVD) technique was used to reduce the multiple input multiple output (MIMO) coupled system to a set of the most relevant control system. A linear-quadratic-integral (LQI) controller was then designed for the reduced order model.
This control approach used is very similar to our rotation control strategy which builds its simplified momentum model from a dynamic model (with some experimentally deduced parameters), uses linearization and model reduction as well (projections on Bessel functions) and builds the corresponding reduced optimal controller (LQG) which contains an integrator, an anti-windup but also an observer (Kalman filter) which enables us to rely on inputs and outputs measurements only while taking feedback actions.
Another method (nonlinear PDE control) consisting of using a backstepping boundary control technique has also been applied to current profile control in \cite{Boyer14}. In this work, a backstepping current profile control algorithm (PI controller) was designed for the DIII-D tokamak. This control design technique provided a systematic method to obtain a boundary feedback law through the transformation of a spatially discretized version of the original system into an asymptotically stable target system with desirable properties.
Motivated by the current profile control, in \cite{Xu11}, the authors consider a 1D parabolic system which is similar to our considered diffusion equation for rotation control. A Proper Orthogonal Decomposition (POD) reduced order model was derived and a reduced order bilinear system was obtained. A convergent successive scheme to compute the solution of a finite-time suboptimal control problem defined for this latter reduced order bilinear system was then designed.
The drawback of this method is that it requires to numerically solve a number of ODEs (Riccati equations) at each iteration of this algorithm whereas our method of linearizing the bilinear term enables us to simply solve the Riccati equation once and reuse the results (gains) through all iterations.
\subsubsection{{Shape control}}
To optimally use the space and to ensure good stabilization in large highly elongated tokamaks, the plasma must be maintained as close as possible to the surrounding walls, thus in addition to a vertical control and a plasma current control, an accurate shape control becomes necessary. Because the coil for ohmic heating, the coil for the vertical and radial fields and the coil for the shaping field are usually the same or partially jointly used, this creates a coupling in the input and output parameters of the control systems and therefore adds more complexity to the problem.
In \cite{Gates05}, plasma shape control using real-time equilibrium reconstruction has been implemented on NSTX. The real-time equilibria provide calculations of the flux at points on the plasma boundary, which are used as input to a shape control algorithm known as isoflux control. The flux at the desired boundary location is compared with a reference flux value, and the difference is used as the basic feedback quantity for the poloidal field coils on NSTX.
More recently, \cite{Kolemen11} gives an overview of the shape control implementations and dynamics studies performed on NSTX, in particular, strike point position and X-point height control.
A PID controller for the strike point was tuned by analyzing the step response of the strike point position to the poloidal coil currents, employing the Ziegler-Nichols method. A system identification of the plasma response to the control inputs was used to build the model. An online automatic relay-feedback PID tuning algorithm was then designed.
\subsection{Kinetic control}
\subsubsection{{Burn control}}
To become an economical alternative energy source, nuclear fusion tokamaks must be capable of operating for extended periods of time in a burning plasma mode characterized by a large value of $Q$, the ratio of fusion power to auxiliary power, in other words, we ideally want more ``power out'' than ``power in''. Achieving and maintaining such conditions requires precise control over the plasma density and temperature.
Modulation of the auxiliary power, modulation of the fueling rate, and controlled injection of impurities are considered as possible actuators for the burn control.
In \cite{Mitarai10}, a PID control law was used to regulate fusion power using the deuterium-tritium fueling rate.
In \cite{Leonov05}, a diagonal multi-input, multi-output linear control scheme for burning plasma kinetics was developed by observing actuator influences during numerical simulations of plasmas.
The approximation of the nonlinear burning plasma model by a linearized one for controller design is a common denominator in many model-based controller designs. The model is linearized, a controller is synthesized using linear techniques, and the resulting design is tested on the original nonlinear model.
These controllers succeed in stabilizing the system in nonlinear simulations against a limited set of perturbations and disturbances.
In \cite{Schuster03}, a nonlinear model involving approximate conservation equations for the energy and particles densities was used to synthesize a nonlinear feedback controller (a nonlinear backstepping algorithm) for burn conditions stabilization. The use of nonlinear control techniques removes the limits imposed by linearization and the resulting controller can accommodate very large perturbations but its implementation and applicability on the PCS are not straightforward and requires several iterations whereas linear controller are standard, simple and more adaptive to other control problems.
\subsubsection{{Rotation control}}
In an operating tokamak, each particle has its own velocity and the net sum of velocities of a particle species is the fluid velocity of that species. This fluid velocity can be separated into components; parallel and perpendicular to the flux surfaces. The velocity perpendicular to a flux surface is called convection, and the velocity parallel to the flux surface is called \textbf{rotation}.
We will consider here the toroidal (parallel) component of the velocity $V_{\phi}$ and its angular frequency $\omega = V_{\phi} / R$ where $R$ is the plasma major radius.
Plasma toroidal rotation and its shear have been recognized as a stabilizing mechanism for magnetohydrodynamic (MHD) instabilities such as the neoclassical tearing mode (NTM) \cite{LaHaye10} where a reduced plasma rotation experimentally destabilizes these NTMs at Lower $\beta$ in the DIII-D, NSTX and JET devices by minimizing the effect of error fields that excite these tearing modes. It also helps prevent the resistive wall modes (RWM) \cite{Garofalo02} where these long-wavelength modes are stabilized by a rapid plasma toroidal rotation.
High plasma rotation can have a significant impact on the plasma confinement time by suppressing energy and particle transport to the walls.
\textbf{Neutral beam injection} (NBI) is the dominant source of momentum and therefore rotation in present-day tokamaks. NBI consists of injecting beams of highly energetic neutral particles into the plasma, heating the plasma through collisions, and naturally transferring momentum.
%The NBI system at DIII-D for example, consists of eight beam-lines, each of which can inject a maximum of $2.5$ MW of power into the plasma, four of them are designed to inject in the same direction as the plasma current, aligned with the magnetic axis, two of them are designed to drive co-current with alignment off-axis, and the last two beams are designed to inject in the opposite direction of the plasma current with on-axis alignment.
The NBI system of NSTX (resp. NSTX-U) considered in this work consists of one set (resp. two sets) of three beams which can each inject a maximum of 2\,MW of power into the plasma. Its configuration is shown in Fig~\ref{nstxu} and its injection is spread throughout the plasma flux surfaces which ensures a rotation drive across the plasma poloidal profile.
Several mechanisms have been developed to control and affect plasma rotation beside the main external source of neutral beams injection as toroidal rotation can be influenced by the intrinsic rotation and the 3D magnetic fields due to MHD instabilities or field errors.
The experimental work done in \cite{Solomon10} for example focuses on investigating mechanisms of driving rotation in fusion plasmas without external momentum input (use of intrinsic rotation and non-resonant magnetic fields). It has been found that the torque from these fields can be enhanced at low rotation, which assists in spinning the plasma from rest, and offers increased resistance against plasma slowing.
In this dissertation, we will focus exclusively on the neutral beams injection as the unique source for driving toroidal rotation. All other sources will not be considered.
Rotation control in tokamaks has been already demonstrated using momentum input from injected neutral beams (NBI) as actuators \cite{Scoville07, Yoshida09}.
\cite{Yoshida09} controls the combined ion temperature profile and toroidal rotation for JT-60U device experimentally using real-time measurement and real-time control system consisting of a tuned PID controller.
In \cite{Scoville07}, simultaneous control of the rotation and stored energy was considered for the DIII-D device. In this work, a model-based control algorithm for simultaneous regulation of plasma rotation and $\beta$ has been developed and casted in linear state space form then combined with a proportional-integral-derivative (PID) transfer function to form a closed loop control algorithm which then has been used by the PCS for experimental testing. Our work will show similarities in controlling both the rotation and the stored energy with \cite{Scoville07} but will use a different methodology (model reduction and an optimal controller). More details in the following section.
\section{Contribution of this dissertation}
The main subject of this thesis is to explore \textbf{toroidal rotation profile control} in tokamaks through its direct application on NSTX and NSTX-U devices.
The methodology used throughout the thesis relies heavily on the mechanical engineering tools that are very commonly used in fluid mechanics and flow control such as model reduction and linear feedback controllers. The idea here is to transpose this knowledge and tools and apply them to control the toroidal rotation of plasma in tokamaks.
The thesis begins with a classical plasma problem as a preliminary case study where this methodology of reduced order modeling and control is applied to demonstrate the effectiveness and applicability of our tools.
We then focus on the main topic which is the control of plasma toroidal rotation in a tokamak, to maintain plasma stability for long-pulse operation. We use experimental measurements from the National Spherical Torus Experiment (NSTX) and two different types of actuation: momentum from injected neutral beams and neoclassical toroidal viscosity generated by three-dimensional applied magnetic fields.
Whether based on the data-driven model obtained for NSTX or pure modeling for NSTX-U, a feedback controller is designed, and predictive simulations using the TRANSP plasma transport code show that the controller drives the plasma rotation to various desired profiles given practical constraints on the actuators and the available rotation measurements. Another application of simultaneously controling the toroidal rotation profile and $\beta_n$ stored energy is shown for NSTX-U as well.
The approach used in this dissertation for NSTX-U is quite similar to the one used in \cite{Scoville07} for DIII-D except that the new and unique aspect of it is the use of non-axisymmetric (three-dimensional) magnetic fields as another actuator in closed-loop feedback control to supplement the neutral beam actuator. Rotation alteration by non-resonant, three-dimensional magnetic fields allows more precise, continuous control of the plasma rotation than NBI, as the momentum delivered by the latter occurs in significantly large, discrete increments.
The modeling and control design of the plasma rotation differ also from \cite{Scoville07}. Starting from a diffusion equation, we first proceed to linearize it, then we apply model reduction on it in order to reduce the dimension of the state of the linear optimal controller (LQG) which is a more sophisticated design as it contains both an observer and an integrator. \cite{Scoville07} does not reduce the model but designs a PID controller tuned through rotation simulations. A model-based controller enables us to tune all the gains of the controller offline and therefore drastically reduce the experimental testing.
Furthermore, using TRANSP to test a reference tracking controller for the rotation profile (and stored energy) has never been done before, and it gives confidence in the control design before experimental testing.
Our strategy for NSTX/NSTX-U rotation control is to obtain practical, low-complexity dynamical models useful for implementing relatively simple controllers for tokamaks that can be implemented in the PCS and easily be applied to experiments. Because our starting equation is bilinear like the one considered in \cite{Xu11}, the linearization strategy enables us to use the standard linear optimal control tools that can adapt to any other type of control.
Finally, in our work, a study of the robustness of the controller to some parameters uncertainties, in particular the perpendicular momentum diffusivity profile $\chi_{\phi}$ and the confinement time $\tau_e$ is performed. This will bring more confidence in the controller's robustness in stability and performance as experimental testing is not currently available for rotation control.
\section{Organization of Part I}
\textbf{Chapter 2:} We will present and explain all the background and the tools used in this thesis for the model reduction and the controller design. This will include the balanced truncation and the Bessel functions decomposition for the reduced-order modeling and the different controller designs such as model-based feedback control design and the observer control design. The methodology consisting of doing a model reduction then building a matching controller will be applied to all our plasma control problems.
\textbf{Chapter 3:} We will present our results on the first application of control on a simple plasma drift waves problem as a preliminary case study where the approach defined in chapter 2 of reduced order modeling and control is applied to demonstrate the effectiveness and applicability of our tools.
\textbf{Chapter 4:} This chapter gets into the main topic which is the control of plasma toroidal rotation in the NSTX and NSTX-U tokamaks.
It uses experimental measurements from the NSTX device but will extend the approach to NSTX-U by relying on numerical models and simulations and eventually will increase the complexity by trying to control the rotation and the stored energy simultaneously.
\textbf{Chapter 5:} We summarize our results and possible future directions for this research.
\chapter{Background: Reduced-order modeling and feedback control for linear time invariant systems}
\label{chapitre2}
\section{Overview and motivation}
From the numerical simulation point of view, research problems in plasmas in tokamaks are very challenging due to complex and highly nonlinear dynamics described by coupled multi-variate differential equations in a multi-parameter operating space.
There are generally no analytical solutions available to describe different phenomena occurring in a plasma due to the complexity of interference of different dynamics at the same time. However there are some model simplifications that enable us to understand some properties better; for instance, considering a plasma as a fluid (gas) enables us to state that the pressure is proportional to the product of density and temperature by the thermodynamic equation $P \propto n T$. Even if these quantities are not homogeneous within plasmas, it gives us the intuition needed to understand that the plasma in the tokamak is hotter and denser at its core thus the pressure is also higher at the core.
Another example is the ideal magnetohydrodynamics (MHD) theory which describes the basic behavior of the plasma as a perfectly conducting fluid without distinguishing the different particles composing this fluid. This approximation is sufficiently accurate to be used as a first approximation in almost every magnetic analysis done for tokamak plasma physics, including studies of instabilities and estimation of the plasma shape and position.
A very general and useful approximation technique is model reduction: a mathematical tool derived from control theory and dynamical systems which allows to extract the key coherent structures of these plasma problems and drastically simplify them. This is what we are going to describe and use in this thesis: reduction and reconstruction methods with their direct application to simulation and control of plasma in tokamaks.
This tool was extensively developed primarily for fluid dynamics problems especially flow control where, as in plasma physics, closed-form analytical solutions are rarely available for engineering applications. Mathematical tools have been developed using ideas from dynamical systems, control theory, and geometric mechanics, in order to extract some key structures and conservation laws while simplifying the original problems.
The idea arose when researchers were focusing on flow control when trying to reduce some undesired properties such as drag or instabilities (vibrations) or to enhance other ones. Major breakthroughs have been made in feedback active control due to the increase of machine capacity in simulations and modeling \cite{ Kim07, Cattafesta08, Choi07, Sipp10} thus the focus shifted to model based feedback control methods but the drawback of these methods was that it was applicable only for limited dimension systems (not bigger than $10^4$) whereas dynamical models in computational fluids dynamics dimensions were typically higher (over $10^5$).
Therefore, having a reduced-order model that captures the main low dimensional dynamical structure (when it exists) that accurately reconstructs the input-output dynamics of the full original model is a popular solution used in flow control to solve the problem of high dimensionality. Based on the reduced model, a feedback controller is designed for the full system. Many applications of this methodology can be found, applied to various problems such as noise reduction in cavity flow \cite{Rowley05} or stabilization of an unstable steady state in \cite{Ahuja10}.
In our case, we consider the plasma in a tokamak as an ideal conducting fluid flowing in a torus and, by analogy, we will apply the same mathematical methods of model reduction and model-based feedback control to plasma problems to show its validity and efficacy.
Reduced order methods enable us to gain computational and experimental time and focus on the main dynamics that will matter for our problems.
The development of engineering tools for automatic regulation and control of a system's behavior never stopped progressing especially in the last several decades.
As this work will only focus specifically on linear time invariant (LTI) control theory, where we assume that plants (i.e., systems to be controlled) and controllers have linear dynamics that do not change with time. Although the LTI assumption is restrictive, the control theory tools specific to it are very standard and easy to apply.
In this chapter, Section~\ref{standard} reviews some standard reduced-order modeling methods for Linear Time Invariant (LTI) systems. Section~\ref{optfed} reviews the standard linear feedback control methods.
These methods are used in Chapter~\ref{chapter7}, Chapter~\ref{rot1&} and Chapter~\ref{rot2&}.
%Subsection~\ref{Besselmeth} highlights the Bessel decomposition method that is used in both Chapters~\ref{rot1&} and \ref{rot2&}.
\section{Standard reduced-order modeling methods}
\label{standard}
In this section, various standard techniques for constructing reduced-order models are briefly reviewed for the model reduction of LTI systems used in this thesis.
\subsection{Projection-based model reduction}
This method involves the projection of a model onto a set of modes and is a widely used approach. It is sometimes called the Petrov-Galerkin projection approach.
We start by defining a stable linear time-invariant state-space system as follows:
\begin{equation}
\label{lti1}
\begin{aligned}
\dot{x}&= A x + B u, \\
y &= C x ,
\end{aligned}
\end{equation}
where $x \in \mathbb{R}^n$ is the state (for instance, the state variables at all grid points of the simulation), $u \in \mathbb{R}^p$ is a vector of inputs (for instance, actuators or disturbances), and $y \in \mathbb{R}^q$ is a vector of outputs (for instance, sensor measurements, or other measurable quantities as linear functions of the state). $A \in \mathbb{R}^{n \times n}$, $B \in \mathbb{R}^{n \times p}$, and $C \in \mathbb{R}^{q \times n}$ are respectively the dynamics, control and sensor matrices.
This system of equations (\ref{lti1}) is asymptotically stable if and only if all the eigenvalues of $A$ are located in the left half plane (i.e., they all have negative real parts). It is neutrally stable if it is stable and any eigenvalue of $A$ is on the imaginary axis (pure imaginary), and it is unstable if any eigenvalue of $A$ is in the right-half plane (positive real part).
The goal of model reduction is to obtain an approximate model that captures the dynamic relationship between the inputs $u$ and the outputs $y$ with the smallest dimension possible, so we can write:
\begin{equation}
\label{linSS3}
\begin{aligned}
\dot{x_r}&= A_r x_r + B_r u, \\
y &= C_r x_r ,
\end{aligned}
\end{equation}
where the reduced state variable $x_r \in \mathbb{R}^r$ and $r \ll n$, $A_r \in \mathbb{R}^{r \times r}$, $B_r \in \mathbb{R}^{r \times p}$, and $C_r \in \mathbb{R}^{q \times r}$ are respectively the corresponding reduced dynamics, control and sensor matrices.
The standard projection approach obtains the reduced order model (\ref{linSS3}) by projecting the model (\ref{lti1}) onto a $r$-dimensional subspace spanned by the columns of a matrix $\Phi_r \in \mathbb{R}^{n \times r}$ along a direction that is orthogonal to a $r$-dimensional subspace spanned by the columns of another matrix $\Psi_r \in \mathbb{R}^{n \times r}$. The bases $\Phi_r$ and $\Psi_r$ are bi-orthogonal:
\begin{equation}
\label{linSS4}
\Psi_r^{\dagger} \Phi_r = I_{r \times r},
\end{equation}
where the dagger $\dagger$ denotes the adjoint operator.
By approximating the original state as $x \approx \Phi_r x_r$ in (\ref{lti1}) and applying the orthogonal property (\ref{linSS4}), we can deduce the reduced-order model of the form:
\begin{equation}
\label{linSS5}
\begin{aligned}
\dot{x_r}&= \left(\Psi_r^{\dagger} A \Phi_r \right) \,x_r + \left( \Psi_r^{\dagger} B \right) \, u, \\
y &= \left( C \Phi_r \right) \,x_r.
\end{aligned}
\end{equation}
This bi-orthogonal projection approach can also be used to reduce nonlinear systems, and if $\Phi_r = \Psi_r $, this method is just an orthogonal Galerkin projection.
We can have various projection-based model reduction methods that use different bases $\Phi_r$ and $\Psi_r$ modes and generate different reduced-order models, but the basic projection idea remains the same for all these techniques.
\subsection{Bessel functions}
\label{Besselmeth}
One of the orthogonal Galerkin projection method that we used in this work for rotation control (Chapters~\ref{rot1&} and~\ref{rot2&}) is choosing the Bessel functions as our projecting basis functions. We then write the approximate state in (\ref{lti1}) as
\begin{equation}
x(\rho,t) \approx \sum_{n=1}^{r} a_n (t) \varphi_n(\rho),
\label{decomp1}
\end{equation}
where $\rho$ and $t$ represent the space and time variables. The basis functions are given by
\begin{equation}
\label{decomp2}
\varphi_n(\rho) = J_0(k_n\rho),\qquad n=1,\ldots,r,
\end{equation}
where $J_0$ denotes the Bessel function of the first kind and $k_n$ denotes the $n$-th root of $J_0$. Furthermore, the basis functions satisfy the orthogonality relation
\begin{equation}
\label{decomp3}
\langle \varphi_n,\varphi_m\rangle = 0,\qquad \text{for $m\ne n$},
\end{equation}
where the inner product is defined by
\begin{equation}
\langle f,g \rangle = \int^1 _0 \rho \, f(\rho) \, g(\rho) \, d\rho.
\end{equation}
Inserting the expansion~\eqref{decomp1} into~\eqref{lti1} then taking the inner product with~$\varphi_m$, and using the orthogonality relation~\eqref{decomp3}, we obtain the following reduced-order system
%
\begin{equation}
\begin{aligned}
\dot a_m &= \sum_{n=1}^r \frac{\langle A \varphi_n, \varphi_m\rangle}{\langle \varphi_m,\varphi_m\rangle} a_n +\frac{ \langle B u,\varphi_m\rangle}{\langle \varphi_m,\varphi_m\rangle} ,\\
y &= C \sum_{n=1}^{r} a_n (t) \varphi_n(\rho)
\end{aligned}
\end{equation}
where $ m=1,\ldots,r$ which is a set of $r$ coupled ordinary differential equations for the coefficients $a_m$. Note that the reduced state becomes a vector of coefficients obtained from the projection of the non reduced state onto the Bessel functions basis.
\subsection{Balanced truncation method}
\label{baltrun}
For systems like (\ref{lti1}), the concepts of controllability and observability can be defined and quantified by a pair of symmetric, positive-semidefinite matrices
\begin{subequations}
\begin{align}
W_c &= \int_0 ^{\infty} e^{At} B B^{\dagger} e^{A^{\dagger}t} dt, \\
W_o &= \int_0 ^{\infty} e^{A^\dagger t} C^\dagger C e^{At} dt,
\end{align}
\label{obscontr}
\end{subequations}
called controllability and observability Gramians. The dagger $\dagger$ denotes the adjoint operator.
The controllability Gramian~$W_c$ provides a measure of the influence of input history on the current state (i.e,\ to what degree each state is excited by inputs), and the observability Gramian~$W_o$ measures the influence of an initial state on future outputs with zero control input (i.e,\ to what degree each state excites future outputs). The larger eigenvalues of the controllability (observability) Gramian correspond to the more controllable (observable) states.
We then define the Hankel norm of a system $G$ (defined by its state space realization (\ref{lti1})) as the maximum ratio, over all input signals, between an output signal norm and the given input signal norm
\begin{equation}
\| G \|_H = \sqrt{ \max \, \text{eig}(W_c W_o) } = \sqrt{ \max \, \text{eig}(W_o W_c) },
\end{equation}
we also define the $j$th Hankel singular value $\sigma_j (G)$ as the square root of the $j$th eigenvalue of $W_o W_c$ (or $W_c W_o$ ) with the ordering
\[
\sigma_1(G) \ge \sigma_2(G) \ge \cdots \ge \sigma_j (G) \ge \cdots.
\]
A balanced truncation involves first a coordinate transformation $T$, called the balancing transformation, that simultaneously diagonalizes the controllability and observability matrices defined by \eqref{obscontr}. That is, under a change of coordinates $x=Tz$, the transformed Gramians become
\begin{equation}
T^{-1} W_c ( T^{-1} )^\dagger = T^{\dagger} W_o T = \Sigma,
\end{equation}
where $\Sigma = \operatorname{diag}(\sigma_1,\ldots,\sigma_n)$ is a diagonal matrix which its diagonal entries are the Hankel singular values ordered so that $\sigma_1 \ge \cdots \ge \sigma_n \ge 0$.
A reduced-order model may then be obtained by truncating the states that are least controllable and observable. That is, if $T = \bigl[ T_1\ \,T_2 \bigr]$, and $x = Tz = T_1 z_1 + T_2 z_2$, then a reduced-order model is obtained by setting $z_2=0$, yielding a model of the form
\begin{equation}
\label{redz2}
\begin{aligned}
\dot{z}_1 &= A_r z_1 + B_r u, \\
y &= C_r z_1 ,
\end{aligned}
\end{equation}
The resulting reduced-order balanced model retains the most controllable and observable states and is therefore suitable for capturing the input-output dynamics of the original system.
Quantitatively the balanced truncation procedure guarantees an a priori analytical upper bound of error between the original system and the reduced-order model. If $G(s)=C(sI-A)^{-1}B$ denotes the transfer function of the system~(\ref{lti1}), and $G_r(s)$ denotes the corresponding transfer function of the approximation~(\ref{redz2}), then
\begin{equation}
\label{ineq1a}
|| G - G_r||_{\infty} < 2 \sum_{k= r+1} ^n \sigma_r.
\end{equation}
In addition, any reduced-order model $G_r$ with $r$ states satisfies
\begin{equation}
\label{ineq2b}
|| G - G_r||_{\infty} > \sigma_{r+1},
\end{equation}
where $ \sigma_{r+1}$ is the first neglected Hankel singular value of $G$. This is a fundamental limitation for any reduced-order model. The two inequalities (\ref{ineq1a}) and~(\ref{ineq2b}) provide error bounds.
More details of this method can be found in \cite{Glover84} and \cite{SandP}.
In Chapter~\ref{chapter7}, we applied the balanced truncation method to the Hasegawa-Wakatani (HW) problem in order to build the reduced order model. However for a certain set of parameters, the HW system exhibited some unstable modes (its dynamics matrix $A$ has some unstable eigenvalues, also called right half plane poles).
Because the balanced truncation method only applies to stable models, we could not apply it directly to our HW problem. An extension for unstable models does exist and consists of first decomposing the model into a stable and unstable submodels.
For instance, we can write a system as a sum of two subsystems
\begin{equation}
G(s) = G_u(s)+ G_s(s),
\end{equation}
where $G_u$ contains all the unstable poles (eigenvalues of $A$ with a non-negative real part), and $G_s$ contains the stable ones (eigenvalues of $A$ with a negative real part).
The balanced truncation can therefore be applied to the stable part of the system $G_s$ to find a reduced order approximation $G_{s \, r}$
which can then be added to the unstable subsystem to obtain an approximation of the full model $G(s)$
\begin{equation}
G(s) = G_u(s)+ G_{s \, r}(s).
\end{equation}
\subsection{Methods for high dimensional systems}
\label{highdim}
The fluid mechanics community has been precursor in developing high dimensional reduced-order modeling methods due to the extremely high state sizes needed for the conversion of partial differential equations to ordinary differential equations: modeling flows of fluids requires state sizes of the order of $10^{\, 6}$ to $10^{\, 9}$. Very high-dimensional model reduction is a current field of study and ongoing research. The overall approach of model reduction can still work for large systems by using the appropriate decomposition tools.
This section briefly describes proper orthogonal decomposition (POD), balanced POD, and the eigensystem realization algorithm (ERA), which became very standard techniques. The methods have \textbf{not been used} in this thesis but are given as examples of what can be done and thus complete the overview of the main existing high dimension model reduction tools. In fact, as seen in some examples of current profile control, \cite{Xu11} uses a POD to build the reduced dynamic model.
\subsubsection{The Proper Orthogonal Decomposition (POD)}
Proper Orthogonal Decomposition or POD also known as the principal component analysis or Karhunen-Loève analysis, has been first used in fluids problems by Lumley \cite{Lumley67, Lumley70}, then has been used widely to study fluid flows, model reduction and control \cite{Sirovich87, Aubry88, Holmes96}.
This method can be applied on both linear and nonlinear systems. It constructs an orthogonal set of modes called POD modes directly from experimental data. The POD model reduction consists of an orthogonal Galerkin projection onto the first $r$ leading modes (modes ranked by their intrinsic energy content).
We start by collecting snapshots from simulation of the full model or experimental data and stacking them by columns in a $X \in \text{R}^{n \times m}$ matrix. $n$ is the dimension of the system, $m$ is the number of snapshot taken and $m \ll n$.
The POD modes are deduced to be the orthonormal eigenvectors of $X^{\dagger}X \in \text{R}^{m \times m}$). Therefore, the orthogonal projection using these modes captures the most energetic modes of the simulation or experiment.
It does not capture the input-output dynamics of the original model as the balanced truncation does, and this is the major drawback as the most important modes for control or measurements are not necessarily the most energetic \cite{Smith05, Ilak08}. Another drawback is that reduced order modeling using POD does not conserve stability for stable linear systems \cite{Smith03}.
\subsubsection{The balanced POD}
Another snapshot-based method was introduced by Rowley \cite{Rowley05} for high dimensional LTI systems with high dimensional inputs and outputs, and it approximates the balanced truncation method developed in subsection \ref{baltrun}.
The BPOD algorithm requires an impulse response of the linear dynamics for each actuator, and an impulse response of the adjoint dynamics for each sensor. (If the number of sensors is large, then a method known as output projection can approximate this process with a smaller number of impulse responses.)
The algorithm computes the SVD of a matrix of inner products between direct and adjoint impulse responses, and uses the decomposition to construct modes that transform the high-dimensional dynamics into a reduced-order approximation of the balanced realization.
The main drawback of this method as well as the balanced truncation method is that it relies on adjoints data that in experimental situation cannot be obtained. If we have a model defined, an adjoint system can be easily deduced, but if we want to rely on experimental snapshots, then these methods are limited and nonapplicable.
\subsubsection{ERA}
The eigensystem realization algorithm (ERA) method is analytically equivalent to the balanced POD \cite{Ma11}, but has the advantage that it does not require adjoint impulse responses. Therefore, unlike balanced POD, it is theoretically possible to use ERA on experimental impulse response data to approximate balanced truncation. Furthermore ERA can be less computationally expensive than balanced POD \cite{Ma11}.
The resulting drawback of the ERA method is thus that ERA does not compute the adjoint modes, which can be useful for observability analysis. Therefore, the choice of method between balanced POD and ERA depends on the modeling goals and availability.
However, both the balanced POD and ERA techniques remain a better choice than the POD method when the dynamics are linear or nearly linear.
\section{Optimal feedback control tools}
\label{optfed}
LTI control branches into two main theories, the classical control theory and modern control theory. The main difference between the two theories is the representation and analytical approaches for dynamical systems: classical control theory relies on frequency-domain representations whereas modern control theory relies on time-domain representations.
These two system representations are related to each other by a simple Laplace transform detailed in what follows.
The general state space realization form of a finite dimensional LTI input-output dynamical system $G$ is given in \eqref{lti1}.
Taking the Laplace transform of system $G$ gives us a new system of equations
%
\begin{equation}
\label{lti2}
\begin{aligned}
s \hat{x}(s) &= A \hat{x}(s) + B \hat{u}(s), \\
\hat{y}(s) &= C \hat{x}(s),
\end{aligned}
\end{equation}
which simplifies into
%
\begin{equation}
\label{lti3}
\hat{y}(s) = G(s) \hat{u}(s)
\end{equation}
where $G(s)$ is the transfer function given by $G(s) = C (sI - A)^{-1} B $.
The early studies about feedback control theory were done in the domain of fluid flows, and have only been using classical control tools, mainly because these controllers are easy to design by hand without high-accuracy models. The tendency and focus have then been shifted towards modern control methods because of its ability to design high performing and robust controllers.
This section focuses on particular modern control concepts that have been used throughout this thesis. In particular, it addresses optimal control methods.
The main idea of optimal control for LTI systems is to design a controller $K$ that will achieve the best control performance for a given set of performance constraints.
This section discusses the linear quadratic regulator (LQR) as the optimal full-state feedback controller and the Kalman filter as the optimal state estimator and reviews the linear quadratic Gaussian (LQG) which combines the LQR and the Kalman filter. We used this type of controller for both the Hasegawa-Wakatani system and the rotation problems.
This control is a standard topic in many textbooks, two good references would be \cite{SandP} and \cite{AandM}. We are providing here a brief overview.
\subsection{ Linear quadratic regulator design}
LQR is a full-state feedback control method that consists of designing an optimal controller gain matrix $K$ such that when the linear state feedback law
\begin{equation}
\label{ltilqr}
u=-K x
\end{equation}
is injected in the input-state part of \eqref{lti1}. The quadratic cost function
\begin{equation}
J = \int_0^{\infty} x^T(t) Q \, x(t) + u^T(t) R \, u(t) \, dt,
\end{equation}
where $Q \in \mathbb{R}^{n \times n}$ is the cost matrix, $Q \ge 0$ (positive semidefinite), $R \in \mathbb{R}^{q \times q}$ is the input cost matrix, $R > 0$ (positive definite), is minimized.
$Q$ and $R$ are thus the weight matrices which determine how the cost function penalizes various components of the state and actuator inputs.
The optimal $K$ is deduced by solving a continuous algebraic Riccati equation. More details can be found in \cite{SandP} and \cite{AandM}.
%for $X \in \mathbb{R}^{n \times n}$ ($X>0$)
%\begin{equation}
%A^T X + X A - X B R^{-1} B^T X + Q = 0,
%\end{equation}
%The LQR optimal gain $K$ is thus deduced as
%\begin{equation}
%K = R^{-1}B^T X.
%\end{equation}
Although the LQR controller has guarantees on the stability and robustness of the closed-loop system, in real-time systems its direct implementation is unfeasible because of the lack of the full knowledge of the entire state $x(t)$, we usually only have point wise measurements. Therefore, it is necessary to design and implement a state observer that estimates the state $x(t)$ based on knowledge of the plant actuators inputs $u(t)$ and sensor outputs $y(t)$.
\begin{figure}[htbp]
\centering
\begin{tikzpicture}[auto, >=latex]
\node (input) [sum] {};
\node (plant) [block, fill=green!20, font=\large, right of=input, node distance=10em] {$\dot{x}(t) = A x(t) + B u(t)$};
\node (output) [coord, right of=plant, node distance=10em] {};
\node (controller) [block, font=\large, below of=plant, yshift=-1em] {$K$};
\draw [connector] (plant) -- node [above] {$x \in \mathbb{R}^n$} (output) |- (controller);
\draw [connector] (controller) -| node [pos=0.99] {$-$} (input);
\draw [connector] (input) -- node [above] {$u \in \mathbb{R}^p$} (plant);
\end{tikzpicture}
\caption{Schematic of a full-state feedback system}
\end{figure}
\subsection{ Observer design}
The observer can be designed using a quadratic estimator known as the Kalman filter. This method is optimal if the errors in representing the state $x(t)$ and the measurements $y(t)$ are stochastic Gaussian processes. Such errors typically arise from inaccuracies in the model, external disturbances, and sensor noise.
A standard linear \textbf{observer} reconstructs a state estimate~$\hat x$, with dynamics given by
\begin{equation}
\dot{\hat{x}}(t) = (A- L C) \hat{x}(t) + B u(t) + L y(t),
\label{obslti}
\end{equation}
where the matrices $A,B$ and~$C$ are the same as those in the system~(\ref{lti1}), and $L$ is a matrix of gains chosen such that the state estimate converges quickly relative to the system's dynamics.
Using our linear model, we design an \textbf{optimal observer} (Kalman filter) to find~$L$.
We introduce two zero-mean Gaussian white noise processes, $w$ the process disturbance and $v$ the sensor noise, with respective covariance matrices $W$ and~$V$, into equations (\ref{lti1}) to obtain
\begin{align}
\dot{x}(t) &= A x(t) + B u(t) + w(t),\\
y(t) &= C x(t) + v(t).
\end{align}
Then the covariance of the error in the state estimate is minimized (assuming the noise models are correct) by setting
\begin{equation}
L = P C^T V^{-1},
\end{equation}
where $P \in \mathbb{R}^{n \times n}$ ($P>0$) is a positive-definite symmetric matrix that solves another algebraic Riccati equation.
%\begin{equation}
%A {P} + P {A}^T - P {C}^T V^{-1} C P + W = 0.
%\end{equation}
The observer generates an estimate of the state from the physics model as represented by the state matrix, the inputs and outputs, and once combined to the feedback controller it forms a linear quadratic Gaussian compensator called LQG.
The LQG controller can stabilize any plant that is detectable and stabilizable, meaning that the actuators can \emph{control} and the sensors can \emph{observe} all unstable eigenmodes of $A$. The controller may not always be robust though, then robust control tools can attempt to improve the robustness of the closed-loop systems.
\begin{figure}[htbp]
\centering
\begin{tikzpicture}[auto, >=latex]
\node (input) [sum] {};
\node (plant) [block, fill=green!20, right of=input, node distance=14em, yshift=1em, label=above:Full linear model] {$\begin{aligned} \dot{x}(t) &= A x(t) + B u(t) + w(t) \\ y(t) &= C x(t) + v(t) \end{aligned}$};
\node (output) [coord, right of=plant, node distance=13.5em] {};
\node (observer) [block, fill=blue!20, below of=plant, xshift=2.5em, yshift=-2.5em, label={[name=label_observer]above:{Kalman Filter}}] {$\dot{\hat x}(t) = A \hat{x}(t) - B u(t) - L (y(t) - C \hat{x}(t))$};
\node (controller) [block, below of=plant, xshift=-11em, yshift=-2.5em, label={[name=label_controller]above:LQR}] {$K$};
\draw [connector] (plant) -- node [above] {$x \in \mathbb{R}^q$} (output) |- (observer);
\draw [connector] (observer) -- node [name=xr, above] {$\hat{x}$} (controller);
\draw [connector] (controller) -| node [pos=0.99] {$-$} (input);
\draw [connector] (input) -- node [name=u, above] {$u \in \mathbb{R}^p$} (input-|plant.west);
\node (noise) [above=0em of u] {$v,w$};
\draw [connector] (noise) -- (noise-|plant.west);
\node (-1) [block, fill=red!20, below of=xr, minimum height=1.5em, minimum width=2.4em] {$-1$};
\node (out) [junction, left=0.3em of controller] {};
\node (in) [coord, right=0.8em of observer, yshift=-0.6em] {};
\draw [connector] (out) |- (-1);
\draw [connector] (-1) -| (in) -- (in-|observer.east);
\node [draw, green!50!black, thick, dashed, fit=(observer) (controller) (-1) (label_observer) (label_controller) (out) (in), label={[green!50!black]below:LQG}, yshift=-0.2em] {};
\end{tikzpicture}
\caption{Schematic of an observer (Kalman filter) connected to a controller (LQR) to form a compensator (LQG)}
\end{figure}
%%%\subsection{Integrator, actuators saturation and anti-windup design}
%%%
%%%We can split the feedback control theory into two main applications, one is disturbance rejection and the other one is reference tracking.
%%%
%%%If the primary goal is tracking a desired reference signal, then we want to minimize the steady state error between the output measurements ($y(t)$) and the corresponding desired values ($y_d$). One way to handle such issue is to use \textbf{integral action}.
%%%
%%%We thus introduce a new state variable~$z$ that is the integral of the error:
%%%\begin{equation}
%%% \dot{z}(t) = y_{d} - y(t) = y_{d} - C x(t).
%%% \label{integral01}
%%%\end{equation}
%%%The overall system (\ref{lti1}) can be then written as
%%%\begin{equation}
%%%\renewcommand\arraystretch{0.8}
%%%\frac{\partial}{\partial t} \begin{bmatrix} x(t) \\ z(t) \end{bmatrix}
%%% = { \begin{bmatrix} A & 0 \\ -C & 0 \end{bmatrix} } \begin{bmatrix} x(t) \\ z(t) \end{bmatrix}
%%% + \begin{bmatrix} B \\ 0 \end{bmatrix} u(t) + \begin{bmatrix} 0 \\ I \end{bmatrix} y_{d}
%%%\label{integral02}
%%%\end{equation}
%%%with a new feedback law designed as
%%%\begin{equation}
%%%u(t) = u_d + K (x_d - x(t)) + K_I \!\!\int (y_d - y(t)) \, dt
%%%\end{equation}
%%%where $K$ is the optimal feedback gain defined through LQR and $K_I$ the optimal integral gain. Both of these gains are obtained by solving a continuous algebraic Riccati equation (CARE) with an extended state that includes the integrator defined in equation (\ref{integral01}).
%%%
%%%A drawback of integral control is that if any actuator value is limited to some range $u\in[u_\text{min},u_\text{max}]$ (as they are in most plasma experiments), then the integrator can accumulate error when the actuator is \emph{saturated}, resulting in poor transient performance, a phenomenon known as \textbf{integrator windup}.
%%%
%%%We can avoid these effects by using a standard \textbf{anti-windup} technique (see, e.g., \cite{AandM, Lewis}), in which one feeds back the difference between the desired value of~$u$ and its actual (possibly saturated) value, as shown in the diagram in Figure~\ref{fig:model13}.
%%%
%%%Figure~\ref{fig:model13} shows the schematic of the overall controller, combining the feedback law~(\ref{ltilqr}) with the observer law~(\ref{obslti}), the integrator~(\ref{integral01}) and the anti-windup approach described above.
%%%
%%%\begin{figure*}
%%%\begin{tikzpicture}[x=0.75cm]
%%% \node (r) {$y_{d}$};
%%% \node[junction, right=1.5 of r] (r in) {};
%%% \node[block, right=2 of r in] (F) {$F$};
%%% \node[block, below=of F] (L) {Observer};
%%% \node[sum, right=of L] (sum feedback) {};
%%% \node[block, right=0.9 of sum feedback] (K) {$K$};
%%% \node[sum, right=5 of F, yshift=3] (sum inputs) {};
%%% \node[junction, right=1 of sum inputs] (before sat) {};
%%% \node[block, right=0.5 of before sat] (sat) {
%%% \begin{tikzpicture}
%%% \draw[very thin] (-.4,0) -- (.4,0) (0,-.25) -- (0,.25);
%%% \draw[very thick] (-.4,-.2) -- (-.2,-.2) -- (.2,.2) -- (.4,.2);
%%% \end{tikzpicture}
%%% };
%%% \node[junction, right=0.5 of sat] (after sat) {};
%%% \node[block, right=1.5 of after sat] (P) {Plant};
%%% \node[junction, right=0.7 of P] (P out) {};
%%% \node[right=0.7 of P out] (y) {$y$};
%%% \node[junction, left=1 of L, yshift=3] (y in) {};
%%% \node[coord, left=1 of L, yshift=-3] (sub y in) {};
%%% \node[coord, label=left:$y$, left=2.4 of y in] (y input) {};
%%% \node[block, above=of F] (Ki) {$K_I$};
%%% \node[sum] (sum lqi) at (Ki -| y in) {};
%%% \node[sum, right=of Ki] (AW out) {};
%%% \node[block, right=0.9 of AW out] (integrator) {$\int$};
%%% \node[sum, above=0.5 of sat] (sum AW) {};
%%% \node[block, above=0.5 of sum AW] (AW) {AW};
%%%
%%% \draw[connector] (r) to (r in) to (F);
%%% \draw[connector] (F.east |- sum inputs) to node [above] {$u_{d}$} (sum inputs);
%%% \draw[connector] (F)[yshift=-12] -| node [right, near end] {$x_{d}$} (sum feedback);
%%% \draw[connector] (sum feedback) to (K);
%%% \draw[connector] (K) -| (sum inputs);
%%% \draw[connector] (sum inputs) to (before sat) to (sat);
%%% \draw[connector] (sat) to (after sat) to ++(down:2.6) -| (sub y in) to (sub y in -| L.west);
%%% \draw[connector] (after sat) to node [above, pos=0.7] {$u$} (P);
%%% \draw[connector] (P) to (P out) to (y);
%%% \draw[connector] (P out) -- ++(down:3.5) -| (y input) to (y in) to (y in -| L.west);
%%% \draw[connector] (L) to node [above] {$\hat x$} node [below, very near end] {$-$} (sum feedback);
%%% \draw[connector] (r in) |- (sum lqi);
%%% \draw[connector] (y in) to node [right, pos=0.95] {$-$} (sum lqi);
%%% \draw[connector] (sum lqi) to (Ki);
%%% \draw[connector] (Ki) to (AW out);
%%% \draw[connector] (AW out) to (integrator);
%%% \draw[connector] (integrator) -| (sum inputs);
%%% \draw[connector] (before sat) |- node [below, very near end] {$-$} (sum AW);
%%% \draw[connector] (after sat) |- (sum AW);
%%% \draw[connector] (sum AW) to (AW);
%%% \draw[connector] (AW) to ++(up:0.6) -| (AW out);
%%%
%%% \draw[green!50!black, ultra thick] (1.4,-2.8) rectangle (15.5,3.1);
%%% \node[green!50!black, anchor=south, font=\large\bfseries] at (8.4,3.1) {Controller};
%%% \end{tikzpicture}
%%%
%%%\caption{Global schematic of the controller that combine a feedforward ($F$), a LQR ($K$), an observer, an integrator $(K_I)$ and an anti-windup (AW).}
%%%\label{fig:model13}
%%%\end{figure*}
%%%
%\section{Robust control tools}
%
%Although the closed-loop systems we obtain using the methods described above are stable in theory, multiple factors could make them unstable in practice, or adversely affect their performance.
%The models on which they are based might be inaccurate, some elements of the dynamics might be missing, or some of the parameters might be off.
%All these inaccuracies are referred as \emph{model uncertainty}.
%Fortunately it is possible to precisely describe this uncertainty, put proper bounds on it, and guarantee that our closed-loop systems stay stable and meet all performance specifications in spite of it. More details can be found in~\cite{SandP}.
%
%In this section, we show how to represent system uncertainty resulting from parameter uncertainty and we give conditions for robust stability and performance.
%We describe robustness analysis for systems with a single input and a single output (SISO) first, then for systems with multiple inputs and multiple outputs (MIMO) which will be our practical case in this disertation.
%
%\subsection{Uncertainty and robustness for SISO systems}
%
%Let $G$ be the model representing the system that we want to control.
%We call it the \textbf{nominal model} and we use it to design a controller.
%Because of various uncertainties, or the desire to keep model complexity low (e.g., linearization, model reduction), it is possible that $G$ is not a very accurate representation of the actual underlying system and that the controller we designed would either be unstable or have poor performance when applied to the actual system.
%The actual system can be seen as a perturbed version $G_p$ of our nominal model $G$.
%However, we do not know exactly where $G_p$ lies compared to $G$ but we can be conservative and make our controller robust to any bounded perturbation with sufficiently large bounds to make sure that the set $\Pi$ of all perturbed models contains $G_p$.
%We say that robust stability (resp. performance) is obtained when all plants in $\Pi$ are stable (resp. performant).
%
%One way to put bounds on the perturbations is to rely on parametric uncertainty.
%It is often the case that instead of a precise value, we only know a range of possible values for some parameters in our systems.
%We can use this structured information to compute analytical bounds on the perturbations of the model.
%
%To greatly simplify the analysis while being strictly conservative about stability and performance, we can use a \emph{norm-bounded} description of uncertainty where $\Pi$ is allowed to contain $\mathcal{H}_\infty$ norm-bounded perturbations of our nominal model $G$.
%There are different ways to integrate this uncertainty into a model but one can get intuition on which one to use by observing the patterns obtained by superimposing the Nyquist plots of many perturbed models.
%Here we will only describe multiplicative uncertainty which is one of the most common.
%The perturbed model can be written as
%\begin{equation} \label{eq:multiplicative}
% G_p = G (1 + \Delta W),
%\end{equation}
%where $W$ is a weighting transfer function, and $\Delta$ is a transfer function satisfying $\|\Delta\|_\infty < 1$.
%A block diagram of the perturbed plant is shown in Figure~\ref{fig:example_multiplicative}.
%
%\begin{figure}[htbp]
% \centering
% \begin{tikzpicture}
% \node (u) {$u$};
% \node[junction, right=0.5 of u] (in) {};
% \node[coordinate, right=1 of in] (below Delta) {};
% \node[block, above=0.5 of below Delta] (Delta) {$\Delta$};
% \node[coordinate, right=1.7 of below Delta] (below W) {};
% \node[block, above=0.5 of below W] (W) {$W$};
% \node[sum, right=0.9 of below W] (sum) {};
% \node[block, right=0.5 of sum] (G) {$G$};
% \node[right=0.5 of G] (y) {$y$};
%
% \draw[connector] (u) to (in) to (sum);
% \draw[connector] (sum) to (G);
% \draw[connector] (in) |- (Delta);
% \draw[connector] (Delta) to (W);
% \draw[connector] (W) -| (sum);
% \draw[connector] (G) to (y);
% \end{tikzpicture}
% \caption{Perturbed model $G_p = G (1 + \Delta W)$}
% \label{fig:example_multiplicative}
%\end{figure}
%
%To find $W$, we first observe that $\Delta W = {G_p}/{G} - 1$.
%Thus we can superimpose the Bode plots of ${G_p}/{G} - 1$ for many perturbed plants and choose a low-order upper bound which will be our $W$.
%
%Next we want to check robust stability. Using the Nyquist criterion we require that for any perturbed model $G_p$, the corresponding loop gain $L_p$ should not encircle the -1 point.
%This translates to the condition
%\begin{equation}
% \| W T \|_\infty < 1,
%\end{equation}
%where $T = \frac{L}{1 + L}$ is the complementary sensitivity function.
%
%Finally we want to check robust performance. We simply require that for any perturbed model $G_p$, a certain performance criterion is satisfied.
%The performance criterion is a set of constraints on the sensitivity function $S = \frac{1}{1 + L}$.
%For instance one can specify a minimum bandwidth, which relates to the speed of the response but also the sensitivity to noise, or an upper bound on peaks values which relates to stability margins.
%These specifications are encoded in a weight $W_P$ (where $P$ stands for \emph{performance}, whereas $p$ stands for \emph{perturbed}).
%The condition translates to
%\begin{equation}
% \max_\omega(| W_P S | + | W T |) < 1.
%\end{equation}
%
%\subsection{Robust Stability and performance analysis for MIMO systems}
%
%The robustness analysis for multiple inputs and multiple outputs (MIMO) systems is very similar but slightly more complex than for SISO systems. In a MIMO system, each input may influence any or all outputs and each output may be influenced by any or all inputs.
%
%A first complication of working with a MIMO system is that all transfer functions are now matrices and as such, multiplication is no longer commutative.
%One consequence is that common SISO transfer functions like the loop gain, the sensitivity, and the complementary sensitivity exist in two variants: input or output depending of the order of multiplication.
%Also, looking at one bode plot per input-output pair can quickly become unwieldy, so instead it is often convenient to look at Bode plots of the singular values of a transfer function matrix.
%Furthermore, upper bounds for robustness and performance are often expressed in terms of the largest singular value $\bar\sigma(H)$ for a transfer function matrix $H$.
%
%% Nominal Stability
%
%Before introducing uncertainty, the first step is to verify that the closed-loop nominal system is stable and define performance specifications.
%For the nominal stability, one can inspect the open loop poles and the Nyquist plot of the system and use the Nyquist criterion to conclude.
%
%% Nominal Performance
%
%Performance is a trade-off between getting a fast response (for instance while tracking the reference signal for a tracking problem) and making sure that measurement noise and other disturbances are not amplified too much by the feedback loop.
%Fortunately, these disturbances enter the system at different points in the loop and sometimes at different frequencies than the reference signal.
%
%We define nominal performance by choosing upper bounds on the transfer functions from these disturbance inputs (sometimes called exogenous outputs) to various points in the feedback loop, typically to the output $y$ and to the input $u$.
%There are two important aspects to consider when choosing the bounds.
%The \emph{static} aspect is given by the \emph{peaks} (maximum values) which are related to the quality of the response.
%The \emph{dynamic} aspect is given by the \emph{bandwidth frequency} which is related to the speed of the response.
%In general, a large bandwidth means better performance since high frequency signals are more easily passed on to the outputs, so the rise time is improved, but it also indicates a system which is sensitive to noise and model uncertainties.
%
%In a MIMO system, typical transfer functions of interest include the input and output sensitivities $S_i$ and $S_o$, the input and output complementary sensitivities $T_i$ and $T_o$, and combinations of the sensitivities with the system $G$ and the controller $K$.
%
%% Uncertainty modeling
%
%As in the previous section, we introduce parametric uncertainty which we translate into multiplicative uncertainty.
%In this example, we use right multiplicative uncertainty, as opposed to left (or inverse) multiplicative uncertainty, since the order of multiplications matters.
%The perturbed plant $G_p$ can be written as
%\begin{equation} \label{eq:multiplicative_MIMO}
% G_p = G (I + W_1 \Delta W_2),
%\end{equation}
%where $W_1$ and $W_2$ are weighting transfer function matrices, and $\Delta$ satisfies $\|\Delta\|_\infty < 1$.
%For simplicity, $W_1$ is often taken to be the identity matrix.
%A block diagram of the perturbed plant is shown on Figure~\ref{fig:example_perturbed_MIMO}.
%