Static and Dynamic Optimization (42111)

Static and Dynamic Optimization (42111)

Build. 303b, room 048 Section for Dynamical Systems

Dept. of Applied Mathematics and Computer Science The Technical University of Denmark

Email: nkpo@dtu.dk phone: +45 4525 3356 mobile: +45 2890 3797

2019-11-10 20:18

Lecture 10: Pontryagins principle

1 / 36

Pontryagins Maximum Principle

Outline of lecture

Recap F9 (End Point Constraints) Free Dynamic Optimization (D+C) End Point Constraints

Constrained Control - Decisions Pontryagins Principle (D) Investment planning Pontryagins Principle (C) Orbit injection (II)

Reading guidance (DO: chapter 4)

2 / 36

Dynamic Optimization (D, free)

Find a sequenceu_i,i= 0, ...,N−1which takes the system

xi+1=fi(xi, ui) x0=x₀ from its initial statex₀ along a trajectory such that the performance index

J=φN[xN] +

N−1

X

i=0

Li(xi, ui) is optimized. Define the Hamiltonian function as:

Hi=Li(xi, ui) +λ^T_i+1fi(xi, ui) The Euler-Lagrange equations:

x_i+1=f_i(xi, u_i) λ^T_i = ∂

∂x_iH_i 0 = ∂

∂ui

Hi

Dynamic Optimization (C, free)

Find a functionutt∈[0; T]which takes the system system

˙

x=ft(xt, ut) x₀=x₀ from its initial statex₀ along a trajectory such that the performance index

J=φ_T[xT] + Z T

0

Lt(xt, ut)dt is optimized. Define the Hamilton function as:

H(x, u, λ, t) =Lt(xt, ut) +λ^T_tft(xt, ut) The Euler-Lagrange equations:

˙

x=ft(xt, ut) −λ˙^T = ∂

∂xt

Ht

0 = ∂

∂ut

Ht

3 / 36

Free DDO

with boundary conditions:

x₀=x₀ λ^T_N= ∂

∂x_Nφ_N(xN)

DDO, End points constraints (EPC)

with boundary conditions:

x0=x₀ ψN(xN) = 0 λ^T_N=ν^T ∂

∂xN

ψN(xN)+ ∂

∂xN

φN(xN)

Free CDO

with boundary conditions:

x₀=x₀ λ^T_T= ∂

∂x_Tφ_T(xT)

CDO, End point constraints (EPC)

with boundary conditions:

x0=x₀ ψT(xT) = 0 λ^T_T=ν^T ∂

∂xT

ψT(xT)+ ∂

∂xT

φT(xT)

4 / 36

End point constraints (EPC)

x₀=x₀ ψ_T(xT) = 0 λ^T_T=ν^T ∂

∂x_Tψ_T(xT) + ∂

∂x_Tφ_T(xT)

Simple EPC

x_T=x_T λ^T_T=ν^T+ ∂

∂x_Tφ_T(xT)

Partial simple EPC

xT= x˜T

¯ xT

˜

x_T= ˜x_T ¯x_T is free The boundary conditions becomes:

˜

xT= ˜x_T ˜λ^T_T=ν^T+ ∂

∂x˜T

φT(xT)

¯

xT is free λ¯T= ∂

∂x¯T

φT(xT)

Linear EPC

CxT =r C: p×n matrix The boundary conditions are:

Cx_T =r λ^T_T =ν^TC+ ∂

∂xT

φ_T(xT)

5 / 36

Pontryagins Maximum principle

Constrained decisions:

ui∈ Uⁱ Example:

|ui| ≤¯u Example:

u≤ui≤u¯ Example:

u_i≤ui≤u¯i

Example:

ui∈ {−1, 0, 1}

6 / 36

Pontryagin

Lev Semenovich Pontryagin (3 September 1908 - 3 May 1988) was a Soviet Russian mathematician. He was born in Moscow and lost his eyesight in a primus stove explosion when he was 14.

He made major discoveries in a number of fields of mathematics, including the geometric parts of topology. Later in his career he worked in optimal control theory.

His maximum principle is fundamental to the modern theory of optimization.

Pontryagin was quite a controversial personality.

Source: Wikipedia

7 / 36

Pontryagin (D)

Find a sequenceui,i= 0, ...,N−1where ui∈ Uⁱ which takes the system

x_i+1=fi(xi, ui) x₀=x₀ from its initial statex₀along a trajectory such that the performance index

J=φ[xN] +

N−1

X

i=0

L_i(xi, u_i)

is optimized (minimized or maximized). Defining the Hamiltonian function H_i=L_i(xi, u_i) +λ^T_i+1f_i(xi, u_i) The necessary equations:

x_i+1=fi(xi, ui) λ^T_i = ∂

∂xi

Hi ui=arg min

ui∈Ui

[Hi]

with boundary conditions:

x0=x₀ λ^T_N= ∂

∂xN

φN(xN) If EPC present the last is as given in Chapter 3.

8 / 36

Example: Investment planning

Plan: During a period of time (N) to invest a amount of moneyu_i(limitted to max 600 Dkr) each interval to obtain a specified sum (xN).

Dynamics:

x_i+1= (1 +α)xi+ui x₀= 0 x_N= 10.000Dkr Objective:

M in J J=

N−1

X

i=0

1 2u²_i subject to:

0≤ui≤600Dkr

9 / 36

1 2 3 4 5 6 7 8 9 10 0

200 400 600 800

Input sequence

1 2 3 4 5 6 7 8 9 10 11

0 2000 4000 6000 8000 10000 12000

Saldo

10 / 36

The Hamiltonian function Hi= 1

2u²_i+λ_i+1[axi+ui] a= 1 +α EL (or KKT) conditions:

x_i+1 = ax_i+u_i x₀= 0 x_N= 10000

λi = aλ_i+1 λ_N=ν

ui = arg min

ui∈Ui

(Hi)

λi=νa^N−i

ui=−λ_i+1 for u≤ui≤¯u (−u≥λ_i+1≥ −u)¯ or

ui=max(u, min(¯u,−νa^N−i−1))

For a givenνsolve the state equation with the control inserted.

Ajustingνsuch that EPC is met

xN=x_N= 10000Dkr

11 / 36

Investment planning with economical (linear) cost

What happens (?) if the objective is changed into:

M in J J=

N−1

X

i=0

ui

In that case:

H=ui+λ_i+1(axi+ui) = (1 +λ_i+1)ui+λ_i+1axi

and Pontryagins principle yields:

xi+1 = axi+ui

λi = aλi+1

ui = arg min

ui∈Ui

(Hi)

As previuos we have the costate evolution (νis a constant or a Lagrange multiplier) λ_i=νa^N−i

The optimization gives:

ui=

u (1 +λ_i+1)>0 λ_i+1>−1

¯

u (1 +λ_i+1)<0 λ_i+1<−1

12 / 36

Pontryagin (C)

Find a functionut t∈[0; T]where

ut∈ U^t which takes the system system

˙

x=ft(xt, ut)

from its initial statex₀along trajectories such that the performance index J=φ_T[xT] +

ZT 0

Lt(xt, ut)dt is optimized. Define the Hamilton function as:

Ht(x, u, λ) =Lt(xt, ut) +λ^T_tft(xt, ut) Then the necessary conditions for this problem can be written as:

˙

x=ft(xt, ut) −λ˙^T= ∂

∂xt

Ht ut=arg min

ut∈Ut

[Ht]

with boundary conditions:

x0=x₀ λT= ∂

∂xT

φT(xT) = ∂

∂xφT

or as in Chapter 3 for EPC.

13 / 36

Orbit injection problem II

θ H

a y

u v

z

The problem is to find the specific thrust force with components,a^z_t anda^y_t, satisfying (a^z_t)²+ (a^y_t)²=a²

such that the terminal horizontal velocity,uT, is maximized subject to the dynamics d

dt







 ut

vt

z y









=







 a^z_t a^y_t ut

vt















 u₀ v0

z0

y0









=







 0 0 0 0









and the terminal constraints

vT= 0 yT=H

J=uT ( φT=uT Lt= 0 )

14 / 36

The Hamilton functions (and others) are

Ht=λ^u_ta^z_t+λ^v_ta^y_t +λ^z_tut+λ^y_tvt φ_T =u_T ψ_T= vT

y_T

= 0

H

The costate equation:

−d dt

λ^u_t λ^v_t λ^z_t λ^y_t

=

λ^z_t λ^y_t 0 0 has the boundary conditions

λ^v_T =νv λ^y_T=νy (fixed state, free costate) λ^u_T = 1 λ^z_T = 0 (free state, fixed costate) resulting in

λ^z_t = 0 λ^y_t =νy

λ^u_t = 1 λ^v_t =νv+νy(T−t)

15 / 36

The maximization of a^z_t

a^y_t

=argmax λ^u_ta^z_t+λ^v_ta^y_t +λ^z_tut+λ^y_tvt

subject to

(a^z_t)²+ (a^y_t)²=a² has the solution:

a^z_t a^y_t

= λ^u_t

λ^v_t

a

p(λ^u_t)²+ (λ^v_t)²

16 / 36

The MP problem

M in b^Tu

st. u^Tu≤a² a≥0 has the solution:

u^∗=− a kbkb

Geometric approach

u

b J

a

17 / 36

Analytic approach

JL=b^Tu+λ(u^Tu−a²) KKT:

u^Tu≤a² b^T+ 2λu^T= 0

u=− b

2λ λ²=b^Tb

4a² u=−b a

√b^Tb

18 / 36

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

Orbit injection problem (TDP)

time (t/T) v

u

a^y a^x

19 / 36

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2

Orbit injection problem (TDP)

y in PU

x in PU

20 / 36

% --- function main2

% ---

T=1; % parameters

a=1;

H=0.2;

parm=[-2.4 4.7]’; % Initial guess on parm x0=zeros(4,1); % Initial state variable opt=optimset; % Options for fsolve opt=optimset(opt,’Display’,’iter’);

parm=fsolve(@erf,parm,opt,T,a,x0,H); % Call fsolve for finding parameters [err,time,xt]=erf(parm,T,a,x0,H); % Call erf ones more for getting the

avt=[]; % state trajectories.

for i=1:length(time), t=time(i);

la=[1; parm(1)+parm(2)*(T-t)];

av=la/sqrt(la’*la)*a; % Thrust force as vector

avt=[avt; av’]; % and stored in a matrix

end

% Here goes the plotting commands. (file: ~nkpo/02711/dist3/main2.m) plot(time,[xt(:,1:2) avt]); grid minor; % Plot

21 / 36

% --- function [err,time,xt]=erf(parm,T,a,x0,H)

% ---

% Determine the end point error (err) given the EPC Lagrange multipliers

% in parm (and the constants that specifies the problem).

Tspan=0:T;

[time,xt]=ode45(@tdpc,Tspan,x0,[],parm,T,a);

xT=xt(end,:)’;

err=[xT(2)-0;

xT(3)-H];

% --- function dx=tdpc(t,x,parm,T,a)

% ---

% System model. Determine the (time) derivative of the state vector

% given the time, state (x) and the EPC Lagrange multipliers.

u=x(1); v=x(2); z=x(3); y=x(4);

p1=parm(1); p2=parm(2);

la=[1; p1+p2*(T-t)];

av=la/sqrt(la’*la)*a; % Specific thrust force as a vector

dx=[av; % remember - a vector

u;

v];

22 / 36

Resource Allocation - investmentplanning, production

Free production

Consider a production

˙

xt=αxt x₀=x₀≥0 whereα >0.

Resource Allocation

Let0≤ut≤1be the fraction kept for production (reinvestment).

Then1−utwill be the fraction to be harvested.

The DO problem is:

˙

xt=αutxt x₀=x₀ xt≥0 and

J= Z T

0

(1−ut)xtdt MaximizeJsubject to0≤ut≤1.

Pontryagin

H=Lt+λ^T_tft= (1−ut)xt+λtαutxt

=xt+ (αλt−1)xtut

˙

xt=αutxt x₀=x₀>0

−λ˙t= 1 + (αλt−1)ut λ_T = 0 ut=

( 1 (αλt−1)xt>0 0 (αλt−1)xt<0

sincext≥0:

ut=











1 λt> 1

α (P roduction) 0 λt< 1

α (Harvest)

23 / 36

Resource Allocation

Harvest

Since

λT = 0

there exist an interval[T1;T] (T1< T) where

λt< 1 α Here (in this interval):

ut= 0

˙

xt= 0 xt=x_T λ˙t=−1 λt= (T−t) From this we have (λT₁ =_α¹ =T−T1)

T1=T −1 α

Production

For0≤t < T1

ut= 1

˙

x=αxt x0=x₀ λ˙t=−αλt λ_T1 = 1

α

xt=x₀e^αt x_T1=x₀e^αT1 λt= 1

α e^α(T1−t))

24 / 36

Resource allocation

Solution summary

T1=T −1 α Then:

ut=

1 for 0≤t < T1

0 for T1< t≤T xt=

x0 e^αt for 0≤t≤T1

x₀ e^αT1 for T₁≤t≤T λt=

₁

α e^α(T1−t)) for 0≤t≤T₁ T−t for T1≤t≤T

t

0 0.5 1 1.5 2

xt

0 1 2 3 4 5 6 7

xt

0 0.5 1 1.5 2

λt

0 1 2 3 4 5 6 7

25 / 36

Resource allocation

t

0 0.5 1 1.5 2

xt

0 1 2 3 4 5 6 7

26 / 36

Resource allocation

x_t

0 0.5 1 1.5 2

λt

0 1 2 3 4 5 6 7

27 / 36

Road construction

28 / 36

Road construction

t

-1 0 1 2 3 4 5 6 7

altitue and slope

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2

2.5 The Terraine and its slope

29 / 36

Road construction

Objective: Find road level,xt, such that J=

ZT t=0

1

2(xt−zt)²dt is minimized. Herezt is the level of terrain.

The dynamic is:

˙

xt=u x0=x₀ where

|ut| ≤a

30 / 36

Road construction

Ht=1

2(xt−zt)²+λtut φ(xT) = 0

˙

xt=u x0=x₀

−λ˙t=xt−zt λT= 0

ut=arg min

|ut|≤a

n1

2(xt−zt)²+λtut

o

31 / 36

λt= Z t

0

(zt−xt)dt Notice:λt= 0forxt=zt.

ut=







a forλ <0

? forλ= 0

−a forλ >0

Optimal trajectories are obtained by concatenation of three types of arcs

Regular arcs whereλt>0andut=−a(maximum downhill slope ars).

Regular arcs whereλt<0andut=a(maximum uphill slope ars).

Singular arcs whereλt= 0and where|ut|< acan take any value.

In that intervalλ˙t= 0and thenxt=zt. Sincex˙=uwe haveu= ˙z.

32 / 36

Assume|z˙|< a

t

-1 0 1 2 3 4 5 6 7

altitue and slope

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2

2.5 The Terraine and its slope

then

xt=zt ut= ˙zt λt= 0

33 / 36

Road construction

34 / 36

Road construction

λ > 0

ut=−a xt=zt₁−a(t−t1) λt=

Z t t₁

(zt−xt)dt

t

< t < t

ut=−a xt₂=zt₁−a(t2−t₁) λt₂=

Z t₂ t₁

(zt−xt)dt

λ

= 0

ut= ˙zt

xt=zt

λt= 0

For determination oft1andt2: Z t₂

t1

(zt−xt)dt= 0 zt₂=zt₁−a(t2−t₁)

35 / 36

Reading guidance

DO: Chapter 4

36 / 36