Retrospective Cost Model Refinement and State Estimation for Space ...
Retrospective Cost Model Refinement and State Estimation for Space Weather Modeling and Prediction Dennis S. Bernstein and Aaron Ridley University of Michigan Thanks to: Ankit Goel Supported by Frederica Darema, AFOSR DDDAS Program 1 1
Modeling for DDDAS model ≠ reality All models are wrong Types of model errors
Parameter errors Errors in dynamics----wrong or missing Everything is uncertain to some extent
But some models are useful! ☺
Some are more useful than others
The “accuracy” of a model is meaningful only relative to its intended purpose Model refinement Model + data = better model 1 2
Matrix Decomposition Trick
𝑘𝑘1 𝑞𝑞1
𝑚𝑚1
𝑘𝑘2
𝑞𝑞2
𝑚𝑚2
𝒌𝒌𝟐𝟐 = 𝒌𝒌𝟐𝟐,𝟎𝟎 + 𝑲𝑲 𝒌𝒌𝟐𝟐,𝟎𝟎 is the nominal value 𝑲𝑲 is the parameter error
1 3
Internal and External Signals 𝒌𝒌𝟐𝟐 = 𝒌𝒌𝟐𝟐,𝟎𝟎 + 𝑲𝑲 𝑲𝑲 is the parameter uncertainty
Measured force
𝑘𝑘1
𝑞𝑞1
𝑤𝑤
Measured displacement
This signal IS measured
This force exerted by the unknown stiffness K is NOT measured
𝑣𝑣
𝑚𝑚1
𝑞𝑞2
𝑘𝑘2
𝑦𝑦 𝑣𝑣
𝑚𝑚2
Force NOT measured
Displacement NOT measured
𝑤𝑤
𝑦𝑦0 Feedback Loop 𝑲𝑲
1
Unknown parameter
This signal IS measured
𝑦𝑦
This distance between the masses is NOT measured
4
Parameter Estimation This signal IS measured
𝑤𝑤
This signal is NOT measured
This signal is computed by retrospective optimization
𝑣𝑣
𝑣𝑣�
𝑘𝑘1 𝑞𝑞1
𝒌𝒌𝟐𝟐,𝟎𝟎 = 𝒌𝒌𝟐𝟐 − 𝑲𝑲
𝑘𝑘2,0
𝑚𝑚1
𝑞𝑞2
𝑚𝑚2
Unknown parameter
𝑲𝑲
Inaccessible
�̇ = 𝑨𝑨𝟎𝟎 𝒙𝒙 � + 𝑭𝑭𝒗𝒗 � + 𝑩𝑩𝑩𝑩 𝒙𝒙 �𝟎𝟎 = 𝒒𝒒 �𝟏𝟏 𝒚𝒚 � = 𝑮𝑮 𝒙𝒙 � 𝒚𝒚
� 𝑲𝑲
� is an estimate of 𝑲𝑲 the unknown parameter K
𝑦𝑦
𝑦𝑦�
𝑦𝑦0
This signal IS measured
This signal is NOT measured
This signal is calculated
𝑦𝑦�0
−
This signal is calculated
+ 𝑧𝑧 = 𝑦𝑦0 − 𝑦𝑦�0 This error signal measures the model mismatch
Parameter update 1 5
A Dynamic and Nonlinear Subsystem 𝑘𝑘1
𝑞𝑞1
Reaction Force
𝑘𝑘2
𝑚𝑚1
𝑚𝑚2
𝑞𝑞2
Main Physical System
Friction
𝑚𝑚3
Force
𝑘𝑘3
Feedback Subsystem Dynamic and nonlinear A feedback subsystem need not be a physically separate entity 1 6
Model Refinement = Adaptive Control! 𝒘𝒘
Main Physical System
𝒗𝒗
Unknown Subsystem
Main Physical System Model
� 𝒗𝒗
𝒘𝒘
𝒚𝒚𝟎𝟎 𝒚𝒚
Unknown Subsystem Model
𝒖𝒖𝟎𝟎
�𝟎𝟎 𝒚𝒚 � 𝒚𝒚
Ideal Plant
Ideal Feedback Controller
Physical Plant
𝒛𝒛 = �𝟎𝟎 𝒚𝒚𝟎𝟎 − 𝒚𝒚
𝒖𝒖
Model Update
𝒚𝒚𝟎𝟎 𝒚𝒚
��0𝟎𝟎 𝒚𝒚 𝑦𝑦 � 𝒚𝒚
Adaptive Feedback Controller Controller Update
𝒛𝒛 = �𝟎𝟎 𝒚𝒚𝟎𝟎 − 𝒚𝒚
1 7
Retrospective Cost Adaptive Control Missile Control • Gains are too aggressive and oscillatory
• Variable forgetting factor makes slight changes in gain evolution • Requires manual tuning
Gains adapting
• Continual gain evolution using KF • Q is chosen by trial and error
1
Retrospective Cost Adaptive Control Aircraft Control with Unknown Aileron Jam
90-deg heading change At t = 200 s, the aileron is jammed at 2 deg
Unknown aileron jam
Rudder compensates!!
Gains adapting
1
Retrospective Cost Adaptive Control Active Interior Noise Control Primary content is suppressed by roughly 18 to 23 dBV from 78 to 117 Hz Spillover occurs in other bands
Goal is to use feedback to suppress road and wind noise
These are broadband disturbances
Classical Bode integral constraint implies that suppression across a certain band necessitates amplification (“spillover”) in other bands.
Idea: Shape the response by frequency-weighting the performance signal 1
Retrospective Cost Model Refinement Circuit Experiment
Current and voltage drops are NOT measured
𝑉𝑉𝑖𝑖𝑖𝑖
Series RLC circuit 𝐶𝐶𝑑𝑑
Driving signal is circuit voltage The only measurement is the voltage drop across the resistor
𝑉𝑉𝑜𝑜𝑜𝑜𝑜𝑜
𝐿𝐿
This signal IS measured
These parameters are unknown and inaccessible
The inductance and capacitance are assumed to be unknown
1
𝒙𝒙 𝑳𝑳𝒙𝒙̈ + 𝑹𝑹𝒙𝒙̇ + = 𝑽𝑽̇𝒊𝒊𝒊𝒊 𝑪𝑪𝒅𝒅 𝑽𝑽𝒐𝒐𝒐𝒐𝒐𝒐 = 𝑹𝑹𝑹𝑹
Retrospective Cost Model Refinement Estimates of L and C True value
Estimated inductance
RCMR estimate Initial estimate
Estimated capacitance
Initial estimate True value RCMR estimate 1
Retrospective Cost Model Refinement Battery Health Monitoring Objective: Monitor battery health by estimating film growth at the negative electrode using charging measurements
𝑤𝑤 𝑘𝑘
𝑢𝑢 𝑘𝑘
Main Battery System
Unknown Film Growth Subsystem
Main Battery System Model 𝑢𝑢� 𝑘𝑘
𝑦𝑦0 𝑘𝑘
Side Reaction Intercalation Current JS
Resistive Film δfilm
𝑦𝑦 𝑘𝑘
𝑦𝑦�0 𝑘𝑘
Adaptive Film Growth Subsystem Model 1𝑦𝑦 �
𝑘𝑘
𝑧𝑧 𝑘𝑘
Retrospective Optimization
= 𝑦𝑦0 − 𝑦𝑦�0
13
Space Weather Modeling
2010
Problem
2012
2013
Unknown changes to the atmospheric density degrade the accuracy of GPS and impede the ability to track space objects
Goals
Use input reconstruction to estimate atmospheric drivers that determine the evolution of the ionosphere-thermosphere Use model refinement to improve the accuracy of atmospheric models Achieve more accurate data assimilation Effects on Earth
Space weather affects the terrestrial environment
Space weather disturbances interfere with satellite and radio communications and operations
Extreme space weather events can knock out the power grid, melt electronics, damage satellites, and disrupt polar air routes
Orbit Determination
Space Debris
1
Orbital prediction error is principally caused by problems in estimating atmospheric drag
Predicting atmospheric drag requires prediction of the atmospheric density and understanding ion-neutral interactions
Measurements in the upper atmosphere are primarily space-based
Monitoring Space Weather
Satellites Solar Missions Magnetospheric Missions Atmospheric Missions
Ground-Based Observatories Ionospheric characteristics and disturbances Atmospheric winds Solar, magnetic, and current indices
Launched 15 July’00 Re-entered 12 Sept’10
Monitoring and Data Centers NOAA Space Weather Prediction Center Heliophysics Events Knowledge base Dominion Observatory in Penticton, British Columbia, Canada Launched 17 March’02 In decaying orbit till 2016 1 Figures: NASA
GITM Global Ionosphere Thermosphere Model Vertical Ion Dynamics
Horizontal Ion Dynamics
Vertical Neutral Dynamics Horizontal Neutral Dynamics
1
GITM Is Fully Parallelized 2x2 Grid of the Earth (4 blocks) Resolution is specified by the number of blocks covering the Earth Each block has 4050 cells 9 longs X 9 lats X 50 alts Each cell as 28 states 7 neutrals, 8 ions, 3 temperatures, 7 neutral velocities, 3 ion velocities Each block has 113,400 states
1
2
3
4
5
6
7
8
9
2 3 4 5
4 0 5 0
6
C
E
L
L
S
1 Block
7 8 9
Typical grids are (longitude blocks)x(latitude blocks): 2x2 (4 processors) Testing purposes, 10° × 20° 453,600 states 8x8 (64 processors) Low resolution physical runs, 5° × 2.5° 7,257,600 states 8x12 (96 processors) High resolution physical runs, 5° × 1.67° 10,886,400 states
1 Block
1 Block
1 17
Estimate Photoelectron Heating Efficiency Vertical Neutral Dynamics
𝑄𝑄EUV = 𝜖𝜖EUV 𝑞𝑞EUV + 𝜖𝜖PHE 𝑞𝑞PHE Photoelectron Heating Efficiency
• 𝑞𝑞EUV is the solar extreme ultraviolet (EUV) heating • 𝑞𝑞PHE is the photoelectron heating---the heat released by secondary photoelectrons 1
Estimate Photoelectron Heating Efficiency Using Artificial Data
Sun
𝑤𝑤
This signal is the measured F10.7
This signal is the retrospectively optimized PHE
𝑣𝑣�
NOTE: This is a purely simulation test of RCMR
Faux Earth
𝑣𝑣
Modeled Modeled PHE PHE
𝑦𝑦
𝑦𝑦�0
This signal is the artificial neutral density along the faux CHAMP orbit
𝑦𝑦0
This signal is the calculated neutral density along the faux CHAMP orbit
Faux CHAMP
−+
𝑧𝑧 = 𝑦𝑦0 − 𝑦𝑦�0 This error signal measures the model mismatch
Estimated Estimated PHE PHE Retrospective Optimization
PHE is an unknown parameter in GITM To estimate PHE, we compute neutral density along CHAMP’s orbit at the fixed altitude of 400 km, and RCMR uses this artificial data As a quality metric for the state estimates,1 we compare estimates of the neutral density along GRACE’s orbit at a constant altitude of 400 km 19
Estimate Photoelectron Heating Efficiency Using Artificial Data
The 90-minute average of the estimated neutral density converges to the artificial neutral density along the CHAMP and GRACE orbits at 400 km altitude
90-minute averaged 𝝆𝝆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆 along the CHAMP orbit
Faux CHAMP neutral density IS assimilated
90 minute averaged 𝝆𝝆𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚 along the CHAMP orbit
Faux CHAMP
90-minute averaged 𝝆𝝆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆 along the GRACE orbit
Faux GRACE neutral density is NOT assimilated
90-minute averaged 𝝆𝝆𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚 along the GRACE orbit
Faux GRACE
1
Faux GRACE neutral density is a performance metric 20
Estimate Photoelectron Heating Efficiency Using Artificial Data
The RCMR estimate of PHE converges to the artificial (modeled) value of PHE
Initial PHE guess Initial PHE guess Modeled PHE
Estimated PHE Modeled PHE Initial Guess Estimated PHE 1 21
Estimate Photoelectron Heating Efficiency Using Artificial Data PHE convergence yields convergent state estimates
1 22
Estimate Photoelectron Heating Efficiency Using Artificial Data
We study the robustness of RCMR to the choice of H The estimate converges to the true value of PHE for a wide range of H
1 23
Estimate Photoelectron Heating Efficiency Using Real Satellite Data
Real Sun
𝑤𝑤
𝑣𝑣
This signal is the measured F10.7
This signal is the retrospectively optimized PHE
𝑣𝑣�
Real Earth Real PHE
𝑦𝑦
𝑦𝑦�0
This signal is the neutral density from real CHAMP
𝑦𝑦0 This signal is the calculated neutral density along the CHAMP orbit
Real CHAMP
−+
𝑧𝑧 = 𝑦𝑦0 − 𝑦𝑦�0 This error signal measures the model mismatch
Estimated Estimated PHE PHE Retrospective Optimization
We estimate PHE using neutral density measurements from the real CHAMP satellite We assimilate real CHAMP satellite data from 2002-11-24 to 2002-12-06 Neutral density measurements from GRACE are used as a quality metric But these data are NOT assimilated 1 Real GRACE
24
Estimate Photoelectron Heating Efficiency Using Real Satellite Data RCMR minimizes the error z in the CHAMP neutral density Real CHAMP neutral density data
GITM+RCMR estimate of neutral density along CHAMP orbit
Real CHAMP
Real CHAMP neutral density IS assimilated 1 25
Estimate Photoelectron Heating Efficiency Using Real Satellite Data
RCMR determines the value of PHE that minimizes the error z in the neutral density estimate
Initial PHE guess
Estimated PHE 1 26
Estimate Photoelectron Heating Efficiency Using Real Satellite Data RCAISE also corrects the neutral density at GRACE’s location GITM+RCAISE estimate of neutral density along GRACE orbit
Real GRACE
GRACE neutral density data is NOT assimilated This data is used only as a quality metric
Δ𝜌𝜌
Real GRACE neutral density data
This error may be due to what scientists believe is a calibration error in the GRACE data
Δ𝜌𝜌 𝑘𝑘 = 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 𝑧𝑧 1: 𝑘𝑘 1
27
Estimate Eddy Diffusion Coefficient in GITM using Total Electron Content
Total electron content (or TEC) is an important descriptive quantity for the ionosphere.
TEC is the total number of electrons integrated between two points, along a tube of one meter squared cross section. Units are 1 TECU=𝟏𝟏𝟏𝟏𝟏𝟏𝟏𝟏 𝒎𝒎−𝟐𝟐
CORS = Continuously Operating Reference Stations. GPS/Met = Ground-Based GPS Meteorology. RTIGS = Real Time International GNSS Service. GNSS = Global Navigation Satellite Systems.
1 TEC plot for the continental USA, made on 11/24/2013
Estimate Eddy Diffusion Coefficient in GITM using Total Electron Content
Solar Flux and tides
𝒘𝒘
𝑻𝑻𝑻𝑻𝑪𝑪𝒎𝒎𝒎𝒎𝒎𝒎𝒎𝒎
GITM with EDC = 1750
This signal is calculated
This signal is computed by retrospective optimization
𝐸𝐸𝐸𝐸𝐶𝐶𝑒𝑒𝑒𝑒𝑒𝑒
GITM with EDC updated every 𝚫𝚫𝒕𝒕 minutes
𝐸𝐸𝐸𝐸𝐶𝐶𝑒𝑒𝑒𝑒𝑒𝑒 = Φ 𝜃𝜃 1
𝑇𝑇𝑇𝑇𝐶𝐶𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐
−
+
𝒛𝒛 = 𝑻𝑻𝑻𝑻𝑪𝑪𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄 − 𝑻𝑻𝑻𝑻𝑪𝑪𝒎𝒎𝒎𝒎𝒎𝒎𝒎𝒎
Estimate Eddy Diffusion Coefficient in GITM using Total Electron Content
EDC=1750 is fixed in the simulated measurement case
EDC estimate computed by RCMR
Initial EDC guess
1
Estimate Eddy Diffusion Coefficient in GITM using Total Electron Content Computed by GITM with RCMR updated EDC
Measurement from simulation with EDC=1750
TEC is measured at a fictional TEC station near Somalia
1
RCMR Block Diagram
𝒘𝒘
𝑻𝑻𝑻𝑻𝑻𝑻𝑻𝑻 𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔
Main System
𝒖𝒖
Unknown Subsystem 𝑮𝑮𝐬𝐬
𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴 𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔
Main system model
� 𝒖𝒖
Subsystem � 𝒔𝒔 Model 𝑮𝑮 1
𝒚𝒚𝟎𝟎 𝒚𝒚
�𝟎𝟎 𝒚𝒚
+
−
� 𝒚𝒚 𝒛𝒛 32
Representational Model with RCMR 𝒘𝒘
True Physical Process 𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴 𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔
Modeled Process
� 𝒖𝒖
Adaptive Subsystem
𝒚𝒚𝟎𝟎 �𝟎𝟎 𝒚𝒚
+
−
� 𝒚𝒚 𝒛𝒛
RCMR can also be used to fit dynamic models to represent unmodeled physics or include subgrid scale effects. 1 33
Representational Model with RCMR Direct numerical simulation
Navier–Stokes equations are numerically solved without a turbulence model. ⇒ all spatial and temporal scales of the turbulence must be resolved. computationally expensive and currently prohibitive for practical problems.
Large eddy simulation (LES)
LES is a mathematical model for turbulence used in CFD. Kolmogorov's (1941) theory of self similarity ⇒ large eddies of the flow are geometry dependent, while smaller scales are more universal. Hence, large eddies can be explicitly solved in a calculation and small eddies are implicitly are accounted for by using a subgrid-scale model (SGS model). Subgrid-scale models
Smagorinsky model (Smagorinsky, 1963) Algebraic Dynamic model (Germano, et. al., 1991) Dynamic Global-Coefficient model (You & Moin, 2007) Localized Dynamic model (Kim & Menon, 1993) 1 WALE (Wall-Adapting Local Eddy-viscosity) model (Nicoud and Ducros, 1999)
34
Representational Model with RCMR 𝒘𝒘
� 𝒖𝒖
DNS or Experiment
𝒚𝒚𝟎𝟎
Large Eddy Simulation
�𝟎𝟎 𝒚𝒚
Subgrid-scale Model
� 𝒚𝒚
+
−
𝒛𝒛
RCMR can also be used to fit dynamic models to represent unmodeled physics or include subgrid scale effects. 1 35
The Burgers Equation
Research with Karthik Duraisamy The Burgers equation is a fundamental PDE in applied mathematics 𝜕𝜕𝜕𝜕 𝜕𝜕 𝑢𝑢2 𝜕𝜕 𝜕𝜕𝜕𝜕 + = 𝜈𝜈 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 2 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕
It is a simplified Navier-Stokes’ momentum equation 𝜕𝜕(𝜌𝜌𝜌𝜌) 𝜕𝜕(𝜌𝜌𝜌𝜌) 𝜕𝜕(𝜌𝜌𝜌𝜌) 𝜕𝜕(𝜌𝜌𝜌𝜌) 𝜕𝜕𝑝𝑝 𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝑣𝑣 𝜕𝜕𝑤𝑤 + 𝑢𝑢 + 𝑣𝑣 + 𝑤𝑤 =− + 𝜈𝜈 + + + 𝜌𝜌𝑔𝑔𝑥𝑥 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝑦𝑦 𝜕𝜕𝑧𝑧
Setting 𝑣𝑣 = 𝑤𝑤 = 𝑝𝑝 = 𝑔𝑔 = 0, 𝜌𝜌 = 1 yields the 1D Burgers equation.
It is a nonlinear diffusive equation, and hence, can qualitatively describe shock waves and viscous diffusion. We apply RCMR to estimate the viscosity in the Burgers equation using the measurement of the flow variable at one point in the domain (grid point 98) 1 Reflects the physical case of limited real measurements 36
Estimate Viscosity in the Burgers Equation
The discretized Burgers equation is 𝑢𝑢𝑗𝑗 𝑘𝑘 + 1 = 𝑢𝑢𝑗𝑗 𝑘𝑘 −
Δ𝑡𝑡 1.5𝑢𝑢𝑗𝑗 𝑘𝑘 2Δ𝑥𝑥
2
− 2𝑢𝑢𝑗𝑗−1 𝑘𝑘
2
+ 0.5𝑢𝑢𝑗𝑗−2 𝑘𝑘
𝑈𝑈 𝑘𝑘 + 1 = 𝐹𝐹 𝑈𝑈 𝑘𝑘
𝑌𝑌0 𝑘𝑘 = 𝑢𝑢98 (𝑘𝑘) 𝑌𝑌 𝑘𝑘 = 𝐻𝐻 𝑈𝑈 𝑘𝑘
+ 𝝂𝝂
Δ𝑡𝑡 𝑢𝑢 (𝑘𝑘) − 2𝑢𝑢𝑗𝑗 (𝑘𝑘) + 𝑢𝑢𝑗𝑗−1 (𝑘𝑘) Δ𝑥𝑥 2 𝑗𝑗+1
T
Defining 𝑈𝑈 𝑘𝑘 = 𝑢𝑢1 𝑘𝑘 𝑢𝑢2 𝑘𝑘 … 𝑢𝑢N 𝑘𝑘
2
+ 𝑊𝑊 𝑘𝑘
𝑊𝑊 𝑘𝑘 = 𝝂𝝂𝑌𝑌 𝑘𝑘
1 37
Estimate Viscosity in Burgers Equation 𝑢𝑢 𝑥𝑥, 0 = 3 + sin(2𝜋𝜋𝜋𝜋) + sin(4𝜋𝜋𝜋𝜋 + 3) + sin(14𝜋𝜋𝜋𝜋 + 5)
Measurement Location 𝑢𝑢98 𝑘𝑘 𝑢𝑢�98 (𝑘𝑘) • • • •
Numerical simulation of the Burgers equation. (a) shows the solution 𝑈𝑈(𝑘𝑘) at 𝑘𝑘 = 100, 200, and 300. (b) shows the measurement 𝑌𝑌0 (𝑘𝑘) = 𝑢𝑢98 (𝑘𝑘) . Note that the measurement converges to a finite value that depends on the initial profile 𝑢𝑢(𝑥𝑥, 0) used in the simulation.
• • 1• • •
Estimation of viscosity using RCMR in the Burgers equation. (a) shows the performance 𝑧𝑧(𝑘𝑘) (b) shows the performance 𝑧𝑧(𝑘𝑘) on logarithmic scale. (c) shows the estimate of viscosity 𝜃𝜃. (d) shows the measurement 𝑢𝑢98 (𝑘𝑘) from the true system as well as 𝑢𝑢�98 (𝑘𝑘)computed by the main system model
38
Modified Burgers Equation We modify Burgers equation by adding a representational subgrid-scale stress. 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕
+
𝜕𝜕 𝑢𝑢2 𝜕𝜕𝜕𝜕 2
=
𝜕𝜕 𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕 𝜈𝜈 𝜕𝜕𝜕𝜕
𝟏𝟏 𝝏𝝏𝝏𝝏 − , 𝟐𝟐 𝝏𝝏𝝏𝝏
where 𝜏𝜏 = −2𝒄𝒄
𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕
We estimate 𝒄𝒄 relating subgrid-scale stresses 𝜏𝜏 to resolved strain rate using measurements of the flow variable at one point in the domain (98th grid point). Again reflects the physical case This example is a precursor to constructing representational models of transport phenomenon in large CFD models to capture unmodeled and subgrid-scale features. 1 39
Estimate Viscosity in Modified Burgers Equation 1 𝑢𝑢 𝑥𝑥, 0 = 4 + ∑6𝑖𝑖=1 sin(2𝜋𝜋𝜋𝜋𝜋𝜋 + 𝑖𝑖 2 ) 6
Estimation of viscosity using RCMR in the modified Burgers equation. • (a) shows the performance 𝑧𝑧(𝑘𝑘) • (b) shows the performance 𝑧𝑧(𝑘𝑘) on logarithmic scale. • (c) shows the estimate of 𝜃𝜃. • (d) shows the measurement 𝑢𝑢98 (𝑘𝑘) from the true system as well as 𝑢𝑢�98 (𝑘𝑘) computed by the main system model.
𝒄𝒄
•
1 40
Latest in RCMR – concurrent optimization
RCMR requires minimal modeling information (H values) Since no analytical model is available, this modeling information has been found by numerical testing For EDC estimation this is tedious We seek an efficient technique for concurrent optimization For adaptive control we have developed concurrent optimization Frantisek Sobolic, Ankit Goel, Dennis S. Bernstein, " Retrospective-Cost Adaptive Control Using Concurrent Controller and Target-Model Optimization," submitted to ACC 2016. 1 41
Concurrent Optimization Setup 𝒘𝒘
𝑻𝑻𝑻𝑻𝑻𝑻𝑻𝑻 𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔
𝒙𝒙𝒌𝒌+𝟏𝟏 =
𝟎𝟎. 𝟑𝟑 𝟎𝟎. 𝟒𝟒 𝟎𝟎 𝒙𝒙𝒌𝒌 + 𝒘𝒘 −𝟎𝟎. 𝟏𝟏 𝟎𝟎. 𝟔𝟔 𝟏𝟏 𝒌𝒌 𝒚𝒚𝟎𝟎 𝒌𝒌 = 𝟏𝟏 𝟎𝟎 𝒙𝒙𝒌𝒌
𝒚𝒚𝟎𝟎
• 𝑧𝑧 is minimized by concurrently optimizing 𝐺𝐺𝑓𝑓 and 𝜃𝜃. • No modeling information is used or guessed in the concurrent optimization.
This parameter is uncertain
�𝒌𝒌+𝟏𝟏 = 𝒙𝒙
� 𝒖𝒖
𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴 𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔
𝟎𝟎. 𝟎𝟎 𝟎𝟎. 𝟒𝟒 𝟏𝟏 𝟎𝟎 �𝒌𝒌 + 𝒙𝒙 𝒖𝒖𝒌𝒌 + 𝒘𝒘 −𝟎𝟎. 𝟏𝟏 𝟎𝟎. 𝟔𝟔 𝟎𝟎 𝟏𝟏 𝒌𝒌 �𝟎𝟎 𝒌𝒌 = 𝟏𝟏 𝟎𝟎 𝒙𝒙 �𝒌𝒌 𝒚𝒚 �𝒌𝒌 = 𝟏𝟏 𝟎𝟎 𝒙𝒙 �𝒌𝒌 𝒚𝒚
� = 𝜽𝜽𝒚𝒚 � 𝒖𝒖
1
� 𝒚𝒚
�𝟎𝟎 𝒚𝒚
+
−
𝒛𝒛
42
Concurrent Optimization
• Window Size = 30. • Filter is 5th order FIR.
These are concurrently optimized by RCMR
𝐻𝐻1 𝐻𝐻2 + 2 +⋯ 𝑞𝑞 𝑞𝑞 • 𝑞𝑞 is the shift operator 𝐺𝐺𝑓𝑓 = 𝐻𝐻0 +
1 43
RCMR/RCAISE versus Standard Methods RCAISE does not provide statistical error measures
No estimate of covariance or probability distribution Uses no priors---not Bayesian Uses no ensemble---only a single simulation
Uses only linear least squares techniques
Requires no adjoint code (none is available for GITM)
Computationally inexpensive
Adds minutes to multi-hour ensemble data assimilation But requires estimates of H’s---determined by numerical testing
May be useful as an adjunct to ensemble codes
For model refinement or input estimation in strongly driven systems----systems whose evolution is primarily due to external inputs
1 44
Summary and Future Work Our goal is to make models more accurate (RCMR) and estimate states and inputs (RCAISE)
RCAISE is applicable to strongly driven systems RCMR refines models by identifying inaccessible subsystems
Future research: Use GITM to search for new physics
Heating/cooling models, eddy diffusion, other physics and data Subgrid modeling?
Key RCMR/RCAISE development goal
Fast tuning of H values for large scale nonlinear applications
RCAISE and RCMR are potentially useful for all applications with unknown inputs and significant modeling errors 1 45
Modeling for DDDAS model ≠ reality All models are wrong Types of model errors
Parameter errors Errors in dynamics----wrong or missing Everything is uncertain to some extent
But some models are useful! ☺
Some are more useful than others
The “accuracy” of a model is meaningful only relative to its intended purpose Model refinement Model + data = better model 1 2
Matrix Decomposition Trick
𝑘𝑘1 𝑞𝑞1
𝑚𝑚1
𝑘𝑘2
𝑞𝑞2
𝑚𝑚2
𝒌𝒌𝟐𝟐 = 𝒌𝒌𝟐𝟐,𝟎𝟎 + 𝑲𝑲 𝒌𝒌𝟐𝟐,𝟎𝟎 is the nominal value 𝑲𝑲 is the parameter error
1 3
Internal and External Signals 𝒌𝒌𝟐𝟐 = 𝒌𝒌𝟐𝟐,𝟎𝟎 + 𝑲𝑲 𝑲𝑲 is the parameter uncertainty
Measured force
𝑘𝑘1
𝑞𝑞1
𝑤𝑤
Measured displacement
This signal IS measured
This force exerted by the unknown stiffness K is NOT measured
𝑣𝑣
𝑚𝑚1
𝑞𝑞2
𝑘𝑘2
𝑦𝑦 𝑣𝑣
𝑚𝑚2
Force NOT measured
Displacement NOT measured
𝑤𝑤
𝑦𝑦0 Feedback Loop 𝑲𝑲
1
Unknown parameter
This signal IS measured
𝑦𝑦
This distance between the masses is NOT measured
4
Parameter Estimation This signal IS measured
𝑤𝑤
This signal is NOT measured
This signal is computed by retrospective optimization
𝑣𝑣
𝑣𝑣�
𝑘𝑘1 𝑞𝑞1
𝒌𝒌𝟐𝟐,𝟎𝟎 = 𝒌𝒌𝟐𝟐 − 𝑲𝑲
𝑘𝑘2,0
𝑚𝑚1
𝑞𝑞2
𝑚𝑚2
Unknown parameter
𝑲𝑲
Inaccessible
�̇ = 𝑨𝑨𝟎𝟎 𝒙𝒙 � + 𝑭𝑭𝒗𝒗 � + 𝑩𝑩𝑩𝑩 𝒙𝒙 �𝟎𝟎 = 𝒒𝒒 �𝟏𝟏 𝒚𝒚 � = 𝑮𝑮 𝒙𝒙 � 𝒚𝒚
� 𝑲𝑲
� is an estimate of 𝑲𝑲 the unknown parameter K
𝑦𝑦
𝑦𝑦�
𝑦𝑦0
This signal IS measured
This signal is NOT measured
This signal is calculated
𝑦𝑦�0
−
This signal is calculated
+ 𝑧𝑧 = 𝑦𝑦0 − 𝑦𝑦�0 This error signal measures the model mismatch
Parameter update 1 5
A Dynamic and Nonlinear Subsystem 𝑘𝑘1
𝑞𝑞1
Reaction Force
𝑘𝑘2
𝑚𝑚1
𝑚𝑚2
𝑞𝑞2
Main Physical System
Friction
𝑚𝑚3
Force
𝑘𝑘3
Feedback Subsystem Dynamic and nonlinear A feedback subsystem need not be a physically separate entity 1 6
Model Refinement = Adaptive Control! 𝒘𝒘
Main Physical System
𝒗𝒗
Unknown Subsystem
Main Physical System Model
� 𝒗𝒗
𝒘𝒘
𝒚𝒚𝟎𝟎 𝒚𝒚
Unknown Subsystem Model
𝒖𝒖𝟎𝟎
�𝟎𝟎 𝒚𝒚 � 𝒚𝒚
Ideal Plant
Ideal Feedback Controller
Physical Plant
𝒛𝒛 = �𝟎𝟎 𝒚𝒚𝟎𝟎 − 𝒚𝒚
𝒖𝒖
Model Update
𝒚𝒚𝟎𝟎 𝒚𝒚
��0𝟎𝟎 𝒚𝒚 𝑦𝑦 � 𝒚𝒚
Adaptive Feedback Controller Controller Update
𝒛𝒛 = �𝟎𝟎 𝒚𝒚𝟎𝟎 − 𝒚𝒚
1 7
Retrospective Cost Adaptive Control Missile Control • Gains are too aggressive and oscillatory
• Variable forgetting factor makes slight changes in gain evolution • Requires manual tuning
Gains adapting
• Continual gain evolution using KF • Q is chosen by trial and error
1
Retrospective Cost Adaptive Control Aircraft Control with Unknown Aileron Jam
90-deg heading change At t = 200 s, the aileron is jammed at 2 deg
Unknown aileron jam
Rudder compensates!!
Gains adapting
1
Retrospective Cost Adaptive Control Active Interior Noise Control Primary content is suppressed by roughly 18 to 23 dBV from 78 to 117 Hz Spillover occurs in other bands
Goal is to use feedback to suppress road and wind noise
These are broadband disturbances
Classical Bode integral constraint implies that suppression across a certain band necessitates amplification (“spillover”) in other bands.
Idea: Shape the response by frequency-weighting the performance signal 1
Retrospective Cost Model Refinement Circuit Experiment
Current and voltage drops are NOT measured
𝑉𝑉𝑖𝑖𝑖𝑖
Series RLC circuit 𝐶𝐶𝑑𝑑
Driving signal is circuit voltage The only measurement is the voltage drop across the resistor
𝑉𝑉𝑜𝑜𝑜𝑜𝑜𝑜
𝐿𝐿
This signal IS measured
These parameters are unknown and inaccessible
The inductance and capacitance are assumed to be unknown
1
𝒙𝒙 𝑳𝑳𝒙𝒙̈ + 𝑹𝑹𝒙𝒙̇ + = 𝑽𝑽̇𝒊𝒊𝒊𝒊 𝑪𝑪𝒅𝒅 𝑽𝑽𝒐𝒐𝒐𝒐𝒐𝒐 = 𝑹𝑹𝑹𝑹
Retrospective Cost Model Refinement Estimates of L and C True value
Estimated inductance
RCMR estimate Initial estimate
Estimated capacitance
Initial estimate True value RCMR estimate 1
Retrospective Cost Model Refinement Battery Health Monitoring Objective: Monitor battery health by estimating film growth at the negative electrode using charging measurements
𝑤𝑤 𝑘𝑘
𝑢𝑢 𝑘𝑘
Main Battery System
Unknown Film Growth Subsystem
Main Battery System Model 𝑢𝑢� 𝑘𝑘
𝑦𝑦0 𝑘𝑘
Side Reaction Intercalation Current JS
Resistive Film δfilm
𝑦𝑦 𝑘𝑘
𝑦𝑦�0 𝑘𝑘
Adaptive Film Growth Subsystem Model 1𝑦𝑦 �
𝑘𝑘
𝑧𝑧 𝑘𝑘
Retrospective Optimization
= 𝑦𝑦0 − 𝑦𝑦�0
13
Space Weather Modeling
2010
Problem
2012
2013
Unknown changes to the atmospheric density degrade the accuracy of GPS and impede the ability to track space objects
Goals
Use input reconstruction to estimate atmospheric drivers that determine the evolution of the ionosphere-thermosphere Use model refinement to improve the accuracy of atmospheric models Achieve more accurate data assimilation Effects on Earth
Space weather affects the terrestrial environment
Space weather disturbances interfere with satellite and radio communications and operations
Extreme space weather events can knock out the power grid, melt electronics, damage satellites, and disrupt polar air routes
Orbit Determination
Space Debris
1
Orbital prediction error is principally caused by problems in estimating atmospheric drag
Predicting atmospheric drag requires prediction of the atmospheric density and understanding ion-neutral interactions
Measurements in the upper atmosphere are primarily space-based
Monitoring Space Weather
Satellites Solar Missions Magnetospheric Missions Atmospheric Missions
Ground-Based Observatories Ionospheric characteristics and disturbances Atmospheric winds Solar, magnetic, and current indices
Launched 15 July’00 Re-entered 12 Sept’10
Monitoring and Data Centers NOAA Space Weather Prediction Center Heliophysics Events Knowledge base Dominion Observatory in Penticton, British Columbia, Canada Launched 17 March’02 In decaying orbit till 2016 1 Figures: NASA
GITM Global Ionosphere Thermosphere Model Vertical Ion Dynamics
Horizontal Ion Dynamics
Vertical Neutral Dynamics Horizontal Neutral Dynamics
1
GITM Is Fully Parallelized 2x2 Grid of the Earth (4 blocks) Resolution is specified by the number of blocks covering the Earth Each block has 4050 cells 9 longs X 9 lats X 50 alts Each cell as 28 states 7 neutrals, 8 ions, 3 temperatures, 7 neutral velocities, 3 ion velocities Each block has 113,400 states
1
2
3
4
5
6
7
8
9
2 3 4 5
4 0 5 0
6
C
E
L
L
S
1 Block
7 8 9
Typical grids are (longitude blocks)x(latitude blocks): 2x2 (4 processors) Testing purposes, 10° × 20° 453,600 states 8x8 (64 processors) Low resolution physical runs, 5° × 2.5° 7,257,600 states 8x12 (96 processors) High resolution physical runs, 5° × 1.67° 10,886,400 states
1 Block
1 Block
1 17
Estimate Photoelectron Heating Efficiency Vertical Neutral Dynamics
𝑄𝑄EUV = 𝜖𝜖EUV 𝑞𝑞EUV + 𝜖𝜖PHE 𝑞𝑞PHE Photoelectron Heating Efficiency
• 𝑞𝑞EUV is the solar extreme ultraviolet (EUV) heating • 𝑞𝑞PHE is the photoelectron heating---the heat released by secondary photoelectrons 1
Estimate Photoelectron Heating Efficiency Using Artificial Data
Sun
𝑤𝑤
This signal is the measured F10.7
This signal is the retrospectively optimized PHE
𝑣𝑣�
NOTE: This is a purely simulation test of RCMR
Faux Earth
𝑣𝑣
Modeled Modeled PHE PHE
𝑦𝑦
𝑦𝑦�0
This signal is the artificial neutral density along the faux CHAMP orbit
𝑦𝑦0
This signal is the calculated neutral density along the faux CHAMP orbit
Faux CHAMP
−+
𝑧𝑧 = 𝑦𝑦0 − 𝑦𝑦�0 This error signal measures the model mismatch
Estimated Estimated PHE PHE Retrospective Optimization
PHE is an unknown parameter in GITM To estimate PHE, we compute neutral density along CHAMP’s orbit at the fixed altitude of 400 km, and RCMR uses this artificial data As a quality metric for the state estimates,1 we compare estimates of the neutral density along GRACE’s orbit at a constant altitude of 400 km 19
Estimate Photoelectron Heating Efficiency Using Artificial Data
The 90-minute average of the estimated neutral density converges to the artificial neutral density along the CHAMP and GRACE orbits at 400 km altitude
90-minute averaged 𝝆𝝆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆 along the CHAMP orbit
Faux CHAMP neutral density IS assimilated
90 minute averaged 𝝆𝝆𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚 along the CHAMP orbit
Faux CHAMP
90-minute averaged 𝝆𝝆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆𝒆 along the GRACE orbit
Faux GRACE neutral density is NOT assimilated
90-minute averaged 𝝆𝝆𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚𝐚 along the GRACE orbit
Faux GRACE
1
Faux GRACE neutral density is a performance metric 20
Estimate Photoelectron Heating Efficiency Using Artificial Data
The RCMR estimate of PHE converges to the artificial (modeled) value of PHE
Initial PHE guess Initial PHE guess Modeled PHE
Estimated PHE Modeled PHE Initial Guess Estimated PHE 1 21
Estimate Photoelectron Heating Efficiency Using Artificial Data PHE convergence yields convergent state estimates
1 22
Estimate Photoelectron Heating Efficiency Using Artificial Data
We study the robustness of RCMR to the choice of H The estimate converges to the true value of PHE for a wide range of H
1 23
Estimate Photoelectron Heating Efficiency Using Real Satellite Data
Real Sun
𝑤𝑤
𝑣𝑣
This signal is the measured F10.7
This signal is the retrospectively optimized PHE
𝑣𝑣�
Real Earth Real PHE
𝑦𝑦
𝑦𝑦�0
This signal is the neutral density from real CHAMP
𝑦𝑦0 This signal is the calculated neutral density along the CHAMP orbit
Real CHAMP
−+
𝑧𝑧 = 𝑦𝑦0 − 𝑦𝑦�0 This error signal measures the model mismatch
Estimated Estimated PHE PHE Retrospective Optimization
We estimate PHE using neutral density measurements from the real CHAMP satellite We assimilate real CHAMP satellite data from 2002-11-24 to 2002-12-06 Neutral density measurements from GRACE are used as a quality metric But these data are NOT assimilated 1 Real GRACE
24
Estimate Photoelectron Heating Efficiency Using Real Satellite Data RCMR minimizes the error z in the CHAMP neutral density Real CHAMP neutral density data
GITM+RCMR estimate of neutral density along CHAMP orbit
Real CHAMP
Real CHAMP neutral density IS assimilated 1 25
Estimate Photoelectron Heating Efficiency Using Real Satellite Data
RCMR determines the value of PHE that minimizes the error z in the neutral density estimate
Initial PHE guess
Estimated PHE 1 26
Estimate Photoelectron Heating Efficiency Using Real Satellite Data RCAISE also corrects the neutral density at GRACE’s location GITM+RCAISE estimate of neutral density along GRACE orbit
Real GRACE
GRACE neutral density data is NOT assimilated This data is used only as a quality metric
Δ𝜌𝜌
Real GRACE neutral density data
This error may be due to what scientists believe is a calibration error in the GRACE data
Δ𝜌𝜌 𝑘𝑘 = 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 𝑧𝑧 1: 𝑘𝑘 1
27
Estimate Eddy Diffusion Coefficient in GITM using Total Electron Content
Total electron content (or TEC) is an important descriptive quantity for the ionosphere.
TEC is the total number of electrons integrated between two points, along a tube of one meter squared cross section. Units are 1 TECU=𝟏𝟏𝟏𝟏𝟏𝟏𝟏𝟏 𝒎𝒎−𝟐𝟐
CORS = Continuously Operating Reference Stations. GPS/Met = Ground-Based GPS Meteorology. RTIGS = Real Time International GNSS Service. GNSS = Global Navigation Satellite Systems.
1 TEC plot for the continental USA, made on 11/24/2013
Estimate Eddy Diffusion Coefficient in GITM using Total Electron Content
Solar Flux and tides
𝒘𝒘
𝑻𝑻𝑻𝑻𝑪𝑪𝒎𝒎𝒎𝒎𝒎𝒎𝒎𝒎
GITM with EDC = 1750
This signal is calculated
This signal is computed by retrospective optimization
𝐸𝐸𝐸𝐸𝐶𝐶𝑒𝑒𝑒𝑒𝑒𝑒
GITM with EDC updated every 𝚫𝚫𝒕𝒕 minutes
𝐸𝐸𝐸𝐸𝐶𝐶𝑒𝑒𝑒𝑒𝑒𝑒 = Φ 𝜃𝜃 1
𝑇𝑇𝑇𝑇𝐶𝐶𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐
−
+
𝒛𝒛 = 𝑻𝑻𝑻𝑻𝑪𝑪𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄 − 𝑻𝑻𝑻𝑻𝑪𝑪𝒎𝒎𝒎𝒎𝒎𝒎𝒎𝒎
Estimate Eddy Diffusion Coefficient in GITM using Total Electron Content
EDC=1750 is fixed in the simulated measurement case
EDC estimate computed by RCMR
Initial EDC guess
1
Estimate Eddy Diffusion Coefficient in GITM using Total Electron Content Computed by GITM with RCMR updated EDC
Measurement from simulation with EDC=1750
TEC is measured at a fictional TEC station near Somalia
1
RCMR Block Diagram
𝒘𝒘
𝑻𝑻𝑻𝑻𝑻𝑻𝑻𝑻 𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔
Main System
𝒖𝒖
Unknown Subsystem 𝑮𝑮𝐬𝐬
𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴 𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔
Main system model
� 𝒖𝒖
Subsystem � 𝒔𝒔 Model 𝑮𝑮 1
𝒚𝒚𝟎𝟎 𝒚𝒚
�𝟎𝟎 𝒚𝒚
+
−
� 𝒚𝒚 𝒛𝒛 32
Representational Model with RCMR 𝒘𝒘
True Physical Process 𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴 𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔
Modeled Process
� 𝒖𝒖
Adaptive Subsystem
𝒚𝒚𝟎𝟎 �𝟎𝟎 𝒚𝒚
+
−
� 𝒚𝒚 𝒛𝒛
RCMR can also be used to fit dynamic models to represent unmodeled physics or include subgrid scale effects. 1 33
Representational Model with RCMR Direct numerical simulation
Navier–Stokes equations are numerically solved without a turbulence model. ⇒ all spatial and temporal scales of the turbulence must be resolved. computationally expensive and currently prohibitive for practical problems.
Large eddy simulation (LES)
LES is a mathematical model for turbulence used in CFD. Kolmogorov's (1941) theory of self similarity ⇒ large eddies of the flow are geometry dependent, while smaller scales are more universal. Hence, large eddies can be explicitly solved in a calculation and small eddies are implicitly are accounted for by using a subgrid-scale model (SGS model). Subgrid-scale models
Smagorinsky model (Smagorinsky, 1963) Algebraic Dynamic model (Germano, et. al., 1991) Dynamic Global-Coefficient model (You & Moin, 2007) Localized Dynamic model (Kim & Menon, 1993) 1 WALE (Wall-Adapting Local Eddy-viscosity) model (Nicoud and Ducros, 1999)
34
Representational Model with RCMR 𝒘𝒘
� 𝒖𝒖
DNS or Experiment
𝒚𝒚𝟎𝟎
Large Eddy Simulation
�𝟎𝟎 𝒚𝒚
Subgrid-scale Model
� 𝒚𝒚
+
−
𝒛𝒛
RCMR can also be used to fit dynamic models to represent unmodeled physics or include subgrid scale effects. 1 35
The Burgers Equation
Research with Karthik Duraisamy The Burgers equation is a fundamental PDE in applied mathematics 𝜕𝜕𝜕𝜕 𝜕𝜕 𝑢𝑢2 𝜕𝜕 𝜕𝜕𝜕𝜕 + = 𝜈𝜈 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 2 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕
It is a simplified Navier-Stokes’ momentum equation 𝜕𝜕(𝜌𝜌𝜌𝜌) 𝜕𝜕(𝜌𝜌𝜌𝜌) 𝜕𝜕(𝜌𝜌𝜌𝜌) 𝜕𝜕(𝜌𝜌𝜌𝜌) 𝜕𝜕𝑝𝑝 𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝑣𝑣 𝜕𝜕𝑤𝑤 + 𝑢𝑢 + 𝑣𝑣 + 𝑤𝑤 =− + 𝜈𝜈 + + + 𝜌𝜌𝑔𝑔𝑥𝑥 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝑦𝑦 𝜕𝜕𝑧𝑧
Setting 𝑣𝑣 = 𝑤𝑤 = 𝑝𝑝 = 𝑔𝑔 = 0, 𝜌𝜌 = 1 yields the 1D Burgers equation.
It is a nonlinear diffusive equation, and hence, can qualitatively describe shock waves and viscous diffusion. We apply RCMR to estimate the viscosity in the Burgers equation using the measurement of the flow variable at one point in the domain (grid point 98) 1 Reflects the physical case of limited real measurements 36
Estimate Viscosity in the Burgers Equation
The discretized Burgers equation is 𝑢𝑢𝑗𝑗 𝑘𝑘 + 1 = 𝑢𝑢𝑗𝑗 𝑘𝑘 −
Δ𝑡𝑡 1.5𝑢𝑢𝑗𝑗 𝑘𝑘 2Δ𝑥𝑥
2
− 2𝑢𝑢𝑗𝑗−1 𝑘𝑘
2
+ 0.5𝑢𝑢𝑗𝑗−2 𝑘𝑘
𝑈𝑈 𝑘𝑘 + 1 = 𝐹𝐹 𝑈𝑈 𝑘𝑘
𝑌𝑌0 𝑘𝑘 = 𝑢𝑢98 (𝑘𝑘) 𝑌𝑌 𝑘𝑘 = 𝐻𝐻 𝑈𝑈 𝑘𝑘
+ 𝝂𝝂
Δ𝑡𝑡 𝑢𝑢 (𝑘𝑘) − 2𝑢𝑢𝑗𝑗 (𝑘𝑘) + 𝑢𝑢𝑗𝑗−1 (𝑘𝑘) Δ𝑥𝑥 2 𝑗𝑗+1
T
Defining 𝑈𝑈 𝑘𝑘 = 𝑢𝑢1 𝑘𝑘 𝑢𝑢2 𝑘𝑘 … 𝑢𝑢N 𝑘𝑘
2
+ 𝑊𝑊 𝑘𝑘
𝑊𝑊 𝑘𝑘 = 𝝂𝝂𝑌𝑌 𝑘𝑘
1 37
Estimate Viscosity in Burgers Equation 𝑢𝑢 𝑥𝑥, 0 = 3 + sin(2𝜋𝜋𝜋𝜋) + sin(4𝜋𝜋𝜋𝜋 + 3) + sin(14𝜋𝜋𝜋𝜋 + 5)
Measurement Location 𝑢𝑢98 𝑘𝑘 𝑢𝑢�98 (𝑘𝑘) • • • •
Numerical simulation of the Burgers equation. (a) shows the solution 𝑈𝑈(𝑘𝑘) at 𝑘𝑘 = 100, 200, and 300. (b) shows the measurement 𝑌𝑌0 (𝑘𝑘) = 𝑢𝑢98 (𝑘𝑘) . Note that the measurement converges to a finite value that depends on the initial profile 𝑢𝑢(𝑥𝑥, 0) used in the simulation.
• • 1• • •
Estimation of viscosity using RCMR in the Burgers equation. (a) shows the performance 𝑧𝑧(𝑘𝑘) (b) shows the performance 𝑧𝑧(𝑘𝑘) on logarithmic scale. (c) shows the estimate of viscosity 𝜃𝜃. (d) shows the measurement 𝑢𝑢98 (𝑘𝑘) from the true system as well as 𝑢𝑢�98 (𝑘𝑘)computed by the main system model
38
Modified Burgers Equation We modify Burgers equation by adding a representational subgrid-scale stress. 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕
+
𝜕𝜕 𝑢𝑢2 𝜕𝜕𝜕𝜕 2
=
𝜕𝜕 𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕 𝜈𝜈 𝜕𝜕𝜕𝜕
𝟏𝟏 𝝏𝝏𝝏𝝏 − , 𝟐𝟐 𝝏𝝏𝝏𝝏
where 𝜏𝜏 = −2𝒄𝒄
𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕
We estimate 𝒄𝒄 relating subgrid-scale stresses 𝜏𝜏 to resolved strain rate using measurements of the flow variable at one point in the domain (98th grid point). Again reflects the physical case This example is a precursor to constructing representational models of transport phenomenon in large CFD models to capture unmodeled and subgrid-scale features. 1 39
Estimate Viscosity in Modified Burgers Equation 1 𝑢𝑢 𝑥𝑥, 0 = 4 + ∑6𝑖𝑖=1 sin(2𝜋𝜋𝜋𝜋𝜋𝜋 + 𝑖𝑖 2 ) 6
Estimation of viscosity using RCMR in the modified Burgers equation. • (a) shows the performance 𝑧𝑧(𝑘𝑘) • (b) shows the performance 𝑧𝑧(𝑘𝑘) on logarithmic scale. • (c) shows the estimate of 𝜃𝜃. • (d) shows the measurement 𝑢𝑢98 (𝑘𝑘) from the true system as well as 𝑢𝑢�98 (𝑘𝑘) computed by the main system model.
𝒄𝒄
•
1 40
Latest in RCMR – concurrent optimization
RCMR requires minimal modeling information (H values) Since no analytical model is available, this modeling information has been found by numerical testing For EDC estimation this is tedious We seek an efficient technique for concurrent optimization For adaptive control we have developed concurrent optimization Frantisek Sobolic, Ankit Goel, Dennis S. Bernstein, " Retrospective-Cost Adaptive Control Using Concurrent Controller and Target-Model Optimization," submitted to ACC 2016. 1 41
Concurrent Optimization Setup 𝒘𝒘
𝑻𝑻𝑻𝑻𝑻𝑻𝑻𝑻 𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔
𝒙𝒙𝒌𝒌+𝟏𝟏 =
𝟎𝟎. 𝟑𝟑 𝟎𝟎. 𝟒𝟒 𝟎𝟎 𝒙𝒙𝒌𝒌 + 𝒘𝒘 −𝟎𝟎. 𝟏𝟏 𝟎𝟎. 𝟔𝟔 𝟏𝟏 𝒌𝒌 𝒚𝒚𝟎𝟎 𝒌𝒌 = 𝟏𝟏 𝟎𝟎 𝒙𝒙𝒌𝒌
𝒚𝒚𝟎𝟎
• 𝑧𝑧 is minimized by concurrently optimizing 𝐺𝐺𝑓𝑓 and 𝜃𝜃. • No modeling information is used or guessed in the concurrent optimization.
This parameter is uncertain
�𝒌𝒌+𝟏𝟏 = 𝒙𝒙
� 𝒖𝒖
𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴𝑴 𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔𝒔
𝟎𝟎. 𝟎𝟎 𝟎𝟎. 𝟒𝟒 𝟏𝟏 𝟎𝟎 �𝒌𝒌 + 𝒙𝒙 𝒖𝒖𝒌𝒌 + 𝒘𝒘 −𝟎𝟎. 𝟏𝟏 𝟎𝟎. 𝟔𝟔 𝟎𝟎 𝟏𝟏 𝒌𝒌 �𝟎𝟎 𝒌𝒌 = 𝟏𝟏 𝟎𝟎 𝒙𝒙 �𝒌𝒌 𝒚𝒚 �𝒌𝒌 = 𝟏𝟏 𝟎𝟎 𝒙𝒙 �𝒌𝒌 𝒚𝒚
� = 𝜽𝜽𝒚𝒚 � 𝒖𝒖
1
� 𝒚𝒚
�𝟎𝟎 𝒚𝒚
+
−
𝒛𝒛
42
Concurrent Optimization
• Window Size = 30. • Filter is 5th order FIR.
These are concurrently optimized by RCMR
𝐻𝐻1 𝐻𝐻2 + 2 +⋯ 𝑞𝑞 𝑞𝑞 • 𝑞𝑞 is the shift operator 𝐺𝐺𝑓𝑓 = 𝐻𝐻0 +
1 43
RCMR/RCAISE versus Standard Methods RCAISE does not provide statistical error measures
No estimate of covariance or probability distribution Uses no priors---not Bayesian Uses no ensemble---only a single simulation
Uses only linear least squares techniques
Requires no adjoint code (none is available for GITM)
Computationally inexpensive
Adds minutes to multi-hour ensemble data assimilation But requires estimates of H’s---determined by numerical testing
May be useful as an adjunct to ensemble codes
For model refinement or input estimation in strongly driven systems----systems whose evolution is primarily due to external inputs
1 44
Summary and Future Work Our goal is to make models more accurate (RCMR) and estimate states and inputs (RCAISE)
RCAISE is applicable to strongly driven systems RCMR refines models by identifying inaccessible subsystems
Future research: Use GITM to search for new physics
Heating/cooling models, eddy diffusion, other physics and data Subgrid modeling?
Key RCMR/RCAISE development goal
Fast tuning of H values for large scale nonlinear applications
RCAISE and RCMR are potentially useful for all applications with unknown inputs and significant modeling errors 1 45
