dynobud-1.4.0.0: examples/DaePendulum.hs

{-# OPTIONS_GHC -Wall #-}
{-# Language TypeFamilies #-}
{-# Language FlexibleInstances #-}
{-# Language DeriveFunctor #-}
{-# Language DeriveGeneric #-}
{-# Language DataKinds #-}
{-# Language PolyKinds #-}

module Main where

import GHC.Generics ( Generic, Generic1 )

import Data.Vector ( Vector )

import Accessors

import Dyno.Vectorize
import Dyno.View.View ( View(..), J )
import Dyno.Solvers
import Dyno.Nlp
import Dyno.NlpUtils
import Dyno.Ocp
import Dyno.DirectCollocation.Formulate ( CollProblem(..), makeCollProblem, makeGuess )
import Dyno.DirectCollocation.Types ( CollTraj' )
import Dyno.DirectCollocation.Dynamic ( toMeta )
import Dyno.DirectCollocation.Quadratures ( QuadratureRoots(..) )

import Dynoplot.Callback ( withCallback )

data PendOcp
type instance X PendOcp = PendX
type instance Z PendOcp = PendZ
type instance U PendOcp = PendU
type instance P PendOcp = PendP
type instance R PendOcp = PendR
type instance O PendOcp = PendO
type instance C PendOcp = PendBc
type instance H PendOcp = None
type instance Q PendOcp = None

data PendX a = PendX { pX  :: a
                     , pY  :: a
                     , pVx :: a
                     , pVy :: a
                     , pTorque :: a
                     } deriving (Functor, Generic, Generic1, Show)
data PendZ a = PendZ { pTau :: a}  deriving (Functor, Generic, Generic1, Show)
data PendU a = PendU { pTorqueDot :: a } deriving (Functor, Generic, Generic1, Show)
data PendP a = PendP { pMass :: a } deriving (Functor, Generic, Generic1, Show)
data PendR a = PendR a a a a a a deriving (Functor, Generic, Generic1, Show)
data PendO a = PendO deriving (Functor, Generic, Generic1, Show)
data PendBc a = PendBc (PendX a) (PendX a) deriving (Functor, Generic, Generic1, Show)

instance Vectorize PendX
instance Vectorize PendZ
instance Vectorize PendU
instance Vectorize PendP
instance Vectorize PendR
instance Vectorize PendO
instance Vectorize PendBc

instance Lookup (PendX ())
instance Lookup (PendZ ())
instance Lookup (PendU ())
instance Lookup (PendO ())
instance Lookup (PendP ())

mayer :: a -> PendX a -> PendX a -> None a -> PendP a -> a
mayer tf _ _ _ _ = tf

lagrange :: Floating a => PendX a -> PendZ a -> PendU a -> PendP a -> PendO a -> a -> a -> a
lagrange _ _ u _ _ _ tf = 1e-3*torque'**2 / tf
  where
    PendU torque' = u

r :: Floating a => a
r = 0.3

pendDae :: Floating a => PendX a -> PendX a -> PendZ a -> PendU a -> PendP a -> a -> (PendR a, PendO a)
pendDae (PendX x' y' vx' vy' torque') (PendX x y vx vy torque)
  (PendZ tau) (PendU uTorque') (PendP m) _ = (residual, outputs)
  where
    residual =
      PendR (x' - vx) (y' - vy)
        (m*vx' + x*tau - fx)
        (m*vy' + y*tau - fy)
        (x*vx' + y*vy' + (vx*vx + vy*vy))
        (torque' - uTorque')
    outputs = PendO

    fx =  torque*y
    fy = -torque*x + m*9.8

pendOcp :: OcpPhase' PendOcp
pendOcp = OcpPhase { ocpMayer = mayer
                   , ocpLagrange = lagrange
                   , ocpQuadratures = \_ _ _ _ _ _ _ -> None
                   , ocpDae = pendDae
                   , ocpBc = bc
                   , ocpPathC = pathc
                   , ocpPathCBnds = None
                   , ocpBcBnds = bcBnds
                   , ocpXbnd = xbnd
                   , ocpUbnd = ubnd
                   , ocpZbnd = PendZ (Just (-200), Just 200)
                   , ocpPbnd = PendP (Just 0.3, Just 0.3)
                   , ocpTbnd = (Just 0.1, Just 5)
                   , ocpObjScale      = Nothing
                   , ocpTScale        = Nothing
                   , ocpXScale        = Just pendXScale
                   , ocpZScale        = Just (PendZ 10)
                   , ocpUScale        = Just (PendU 50)
                   , ocpPScale        = Just (PendP 0.3)
                   , ocpResidualScale = Nothing
                   , ocpBcScale       = Just $ PendBc pendXScale pendXScale
                   , ocpPathCScale    = Just None
                   }
pendXScale :: PendX Double
pendXScale = PendX 0.3 0.3 1 1 10

pathc :: Floating a => PendX a -> PendZ a -> PendU a -> PendP a -> PendO a -> a -> None a
pathc _ _ _ _ _ _ = None

xbnd :: PendX Bounds
xbnd = PendX { pX =  (Nothing, Nothing)
             , pY =  (Nothing, Nothing)
             , pVx = (Nothing, Nothing)
             , pVy = (Nothing, Nothing)
             , pTorque = (Just (-30), Just 30)
             }

ubnd :: PendU Bounds
ubnd = PendU (Just (-100), Just 100)

bc :: Floating a => PendX a -> PendX a -> None a -> PendP a -> a -> PendBc a
bc x0 xf _ _ _ = PendBc x0 xf

bcBnds :: PendBc Bounds
bcBnds =
  PendBc
  (PendX
   { pX = (Just 0, Just 0)
   , pY = (Just (-r), Just (-r))
   , pVx = (Just 0, Just 0)
   , pVy = (Just 0, Just 0)
   , pTorque = (Nothing, Nothing)
   })
  (PendX
   { pX = (Nothing, Nothing) -- LICQ
   , pY = (Just r, Just r)
   , pVx = (Just 0, Just 0)
   , pVy = (Nothing, Nothing) -- LICQ
   , pTorque = (Nothing, Nothing)
   })

type NCollStages = 120
type CollDeg = 3

guess :: J (CollTraj' PendOcp NCollStages CollDeg) (Vector Double)
guess = cat $ makeGuess Radau tf guessX guessZ guessU parm
  where
    tf = 1
    guessX t = PendX { pX =   r * sin q
                     , pY = - r * cos q
                     , pVx = r*w*cos q
                     , pVy = r*w*sin q
                     , pTorque = 0
                     }
      where
        q = pi*t/tf
        w = pi/tf
    guessZ _ = PendZ {pTau = 0}
    guessU _ = PendU {pTorqueDot = 0}
    parm = PendP 0.3

solver :: Solver
solver = ipoptSolver { options = [ ("expand", Opt True)
                                 , ("linear_solver", Opt "ma86")
                                 , ("ma86_order", Opt "metis")
                                 ]}

solver2 :: Solver
solver2 = ipoptSolver { options = [("expand", Opt True)] }


main :: IO ()
main = do
  cp  <- makeCollProblem Legendre pendOcp guess
  withCallback $ \send -> do
    let nlp = cpNlp cp
        meta = toMeta (cpMetaProxy cp)
        cb' traj = do
          plotPoints <- cpPlotPoints cp traj
          send (plotPoints, meta)
    _ <- solveNlp solver nlp (Just cb')
--  _ <- solveNlp solver2 nlp Nothing
    return ()